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ABSTRACT 

In the parts of the solar corona and solar wind that experience the fewest Coulomb collisions, the compo- 
nent proton, electron, and heavy ion populations are not in thermal equilibrium with one another. Observed 
differences in temperatures, outflow speeds, and velocity distribution anisotropies are useful constraints on 
proposed explanations for how the plasma is heated and accelerated. This paper presents new predictions of 
the rates of collisionless heating for each particle species, in which the energy input is assumed to come from 
magnetohydrodynamic (MHD) turbulence. We first created an empirical description of the radial evolution of 
Alfven, fast-mode, and slow-mode MHD waves. This model provides the total wave power in each mode as a 
function of distance along an expanding flux tube in the high-speed solar wind. Next we solved a set of cascade 
advection-diffusion equations that give the time-steady wavenumber spectra at each distance. An approximate 
term for nonlinear coupling between the Alfven and fast-mode fluctuations is included. For reasonable choices 
of the parameters, our model contains enough energy transfer from the fast mode to the Alfven mode to excite 
the high-frequency ion cyclotron resonance. This resonance is efficient at heating protons and other ions in the 
direction perpendicular to the background magnetic field, and our model predicts heating rates for these species 
that agree well with both spectroscopic and in situ measurements. Nonetheless, the high-frequency waves com- 
prise only a small part of the total Alfvenic fluctuation spectrum, which remains highly two-dimensional as is 
observed in interplanetary space. 

Subject headings: magnetohydrodynamic s (MHD) — plasmas — solar wind — Sun: corona — turbulence — 
waves 



1. INTRODUCTION 

The energy that heats the solar corona and accelerates the 
solar wind originates in convective motions beneath the Sun's 
surface. However, even after many years of investigation, the 
physical processes that transport a fraction of this energy to 
the corona and convert it into thermal, magnetic, and kinetic 
energy are still not understood. In order to construct and test 
theoretical models, a wide range of measurements of rele- 
vant plasma parameters must be available. In the low-density, 
open-field regions that reach into interplanetary space, the 
number of plasma parameters that need to be measured is 
larger because the plasma becomes collisionless and individ- 
ual particle species (e.g., protons, electrons, and heavy ions) 
can exhibit divergent properties. Such differences in particle 
velocity distributions are valuable probes of kinetic processes 
of heating and acceleration. 

The spectroscopic instruments aboard the Solar and He- 
liospheric Observatory (SOHO) — e.g., the Ultraviolet Coro- 
nagraph Spectrometer (UVCS) and Solar Ultraviolet Mea- 
surements of Emitted Radiation (SUMER) — have measured 
several key collisionless plas ma properties fo r a variety of 
solar wind source regions dKohl et alJ I1995L 1 19971 120061: 
IWilhelm et al.1 11995L 11997th These observations augment 
decades of in situ plasma and field measurements that show 
similar departures fr om thermal equilib r ium in the collision- 
less solar wind (e.g ., lNeugebaudll98l lMarschlll999l 120061: 
Kasper et al. 2008|). In the high-speed solar wind, both coro- 
nal and heliospheric measurements point to the existence of 
preferential ion heating and acceleration, as well as protons 
being hotter than electrons. There are also marked departures 
from Maxwellian velocity distributions for protons and other 
ions, with the temperature measured in directions perpendicu- 
lar to the background magnetic field often exceeding the tem- 



perature parallel to the field (i.e., T± > Tn). 

A large number of different processes have been suggested 
to explain the measured proton and ion properties. Many of 
these processes are related to the dissipation of magnetohy- 
drodynamic (MHD) waves, and many involve multiple steps 
of energy conversion between waves, reconnection structures, 
and other nonlinear plasma features. It was noticed sev- 
eral decades ago that the damping of ion cyclotron resonant 
Alfven waves could naturally give rise to many of the ob- 
serve d plasma propert ies (see reviews by lHollweg & Isenbergl 
2002| iHollwegl l2008l) . The problem in the solar corona, 
though, is how these extremely high-frequency (10 2 -10 4 Hz) 
waves could be generated from pre-existing MHD fluctua- 
tions that appear to have much lower frequencies (< 0.01 Hz). 

One likely source of high-frequency waves and ki- 
netic dissipation is an MHD turbulent cascade. There 
is ample evidence that turbulence provides substantial 
heat input to the plasma in interplanetary space (see 
[ Coleman! \l96§. iGoldstein et all [1991 iTu & Marschl [l99l 
Matth aeus et alJ 120031) . Furthermore, self-consistent mod- 
els of turbulence-driven coronal heating and solar wind ac- 
celeration have begun to succeed in reproducing a wide 
range of observati ons without the need f o r ad hoc free 
parameters (e.g.. ISuzuki & Inutsukal 120061: ICranmer et al. l 
2007 1 : rRappazzo et alj|2008t iBreech et alj|2008t IVerdini et al l 
20ld lBingert & Peter! |20l U Ivan Ballego orTen et al.l l20lTt 
Chandran et al. 1 1201 lh . The general scenario is that convec- 
tion jostles open magnetic flux tubes that are rooted in the 
photosphere and produces Alfven waves that propagate into 
the corona. These waves undergo partial reflection, and the 
resulting "colliding wave packets" drive a turbulent cascade 
which heats the plasma when the eddies reach small enough 
spatial scales. 
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It has been known for many years that Alfvenic turbu- 
lence in a strong magnetic field produces a cascade to small 
scales mai nly in the two-dimensional plane perpendicular t o 
the field (Montgomery & Tumedri98lHShebarin et al.lll983l) . 
and thus is not likely to produce high-frequency ion cyclotron 
waves. In other words, MHD turbulence leads to eddies with 
large perpendicular wavenumbers k± and not large parallel 
wavenumbers ku. Under typical plasma conditions in the 
corona and inner heliosphere, the linear dissipation of high-£^ 
Alfven wa ves would lead to the preferential parallel heating of 
electrons (lLeamonetal.lj l999; Cr anmer & van B allegooiien 
2003; Gary & Borovsky 2008). This apparently disagrees 
with the observational evidence for perpendicular heating of 
positive ions. 

There have been several proposed solutions to the ap- 
parent incompatibility between the predictions of MHD 
turbule nce and existing measurements (see also ICranmerl 
2009a). For example, turbulent fluctuations may be suscep- 
tible to various instabilities that cause ion cyclotron wa ves 
to grow dMarkovskii et alJl2006t IVranjes & Poedtsll2008l) or 
they may induce stochastic perpendic ular motions in ions 
if they reach nonlinear magnitudes dVoitenko & Goossensl 
l200l IWu & ^ng1l2007t IChandranll2010h . Nonetheless, he- 
liospheric measurements have provided several pieces of evi- 
dence for the existence of ion cyclotron resonance that gives 
rise to perpendicular ion heating in the solar wind (e.g., 
Marsch & Tull2001bl: iBourouaine et alJl20Tot iHe et alJl201lt 
Smith et aD,l2012l) . The most direct solution to the problem 
still appears to be for turbulence to transport some fraction of 
the fluctuation energy to high-A;|| cyclotron resonant waves. 

The goal of thi s paper is to investigate the idea proposed by 
IChandrarJ d2005l) for the turbulent generation of ion cyclotron 
waves. In this scenario, nonlinear couplings between Alfven 
waves and other modes such as fast magnetosonic waves pro- 
duce an enhancement in the highly power-law tail of the 
Alfvenic fluctuation spectrum. This is made possible by the 
ability of fast-mode waves to cascade nearly isotropically in 
wavenumber space. Thus, the gradual nonlinear generation 
of ion cyclotron waves may provide enough heat to protons 
and other ions in the corona and inner solar wind (see als o 
iLuo & Melrosell2006t IChandranll2008al: lYoon & Fangll2008l) . 

We note that it is not currently possible to produce a rig- 
orous model that contains a fully self-consistent description 
of MHD wave transport (from the corona to 1 AU), turbulent 
cascade, mode coupling, and dissipation. In order to make 
some progress in trying to understand this complex system, 
we have created models that include a range of simplifying as- 
sumptions. One key approximation is that we divide the mod- 
eling into two separate components: (1) a large-scale model 
of the radial dependence of fluctuation energy densities, and 
(2) a small-scale description of how the "local" fluctuations at 
each radius evolve in wavenumber space and heat the plasma. 
Feedbacks from the latter to the former are not included, and 
we discuss their potential importance in Section [7] 

We model the plasma conditions in a representative mag- 
netic flux tube that is rooted in a polar coronal hole and that 
exhibits a steady-state fast solar wind outflow. In Section|2]we 
describe a model of background plasma conditions and large- 
scale wave transport in this flux tube. We take an empirical ap- 
proach to the solar generation of Alfven, fast, and slow mode 
MHD waves by specifying their amplitudes as free parameters 
at a lower coronal boundary height of 0.01 solar radii (Rq) 
above the photosphere. Section [3] gives a summary of how 
we model the small-scale transport of cascading wave energy 



in wavenumber space, and Section |4] describes our treatment 
of the nonlinear coupling between high-frequency Alfven and 
fast-mode waves. In Section [5] we apply quasilinear kinetic 
theory to predict the net rates of particle heating from the cas- 
cading waves. Section|6]presents a selection of results for the 
collisionless rates of proton, electron, and heavy ion heating. 
Finally, Section|7]concludes this paper with a brief summary 
of our major results, a discussion of some of the wider im- 
plications of this work, and suggestions for future improve- 
ments. 

2. LARGE-SCALE MODEL OF CORONAL HOLE CONDITIONS 

We wish to better understand the global energy budget 
of MHD waves and turbulence from the lower solar corona 
out to the interplanetary medium. The work of this sec- 
tion builds on many earlier models of the radi al evolu- 
tion of Alfven waves in the fast solar wind (e.g., iHolrwegl 
| 1986tlTu & Marschll995HCranmer & van Ballegooijenl2005t 
Chandran & Hollweg l2009ir and extends it to describe the 
likely behavior of fast and slow magnetosonic waves as well. 
Below, we describe an empirical model of how the time- 
steady plasma properties vary with heliocentric distance (Sec- 
tion !2. 11 > as well as a large-scale view of the dispersion, propa- 
gatio n, and di ssipation of linear waves in such a system (Sec- 
tions E2rE!). 

2. 1. Background Time-Steady Plasma 

We model the plasma properties along an open magnetic 
flux tube rooted in a polar coronal hole. At solar mini- 
mum, large unipolar coronal holes are associated with su- 
perradially expanding magnetic fields and the acceleration 
of the high-speed solar wind. Because we only consider 
a field line along the polar axis of symmetry, we do not 
need to include the rotational generation of azimuthal mag- 
netic fields (e.g., the Parker sp iral effect; see IWeber & Davisl 
1 1 967t IPriest & Pneumanll9"74h or other geometri cal effects of 
streamer-like flux tube curvature dLi et al.ll2.01 lb . We do not 
distinguish between dense polar plumes and the more tenu- 
ous interplume regions between them. The radial dependence 
of plasma parameters is described as a function of either the 
heliocentric distance (r) or the height above the solar photo- 
sphere (z= r-R Q ). 

To specify the radial variation of the time-steady mag- 
netic field strength Bo, mass density po, and solar wind 
outflow speed up, we used the empirical description of 
ICranmer & van Ballegooiien (2005). This model combined 
a broad range of observational constraints with a two- 
dimensional magnetostatic model of the expansion of thin 
photospheric flux tubes into a supergranular network canopy. 
At r = 1 AU in this model, the solar wind outflow speed uq is 
78 1 km s" 1 and the proton density n p is 2.56 cm -3 . This model 
also specifies the Alfven speed Va = Bo/(47rpo) 1 ^ 2 , which de- 
creases from a maximum value of 2890 km s" 1 at r = 1.53 Rq 
down to 3 1 km s" 1 at 1 AU. There is a local minimum in Va at 
r ss 1 .02/? Q that is the result of the assumed shape of network 
"funnels" that expand superradially into the corona. 

We also need to know the plasma temperature T in order 
to determine the relative importance of gas pressure versus 
magnetic pressure. Despite observational evidence for dif- 
ferent particle species having different temperatures (and de- 
partures from Maxwellian velocity distributions), we gener- 
ally assume that the majority proton-electron magnetofluid 
is close enough to thermal equilibriu m that strong plasma 
microinstabilities are not excited (e.g.. lGary||199U iMarschl 
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FIG. 1 . — Radial dependence of the Alfven speed (black solid curve), solar 
wind outflow speed (blue dashed curve), one-fluid sound speed (red dotted 
curve), and the angle-averaged, intertial-frame group velocity of fast-mode 
waves (violet dot-dashed curve), all in units of km s . Also shown is the 
dimensionless plasma parameter (green solid curve). 

2006). Thus, we specify a one-fluid temperature T that is 
assumed to be equal to both the proton temperature T p and the 
electron temperature T e , and we assume temperature isotropy 
(Tn Ri 7j_) for both species. 

We used the polar coronal hole model of Cranm er" et al.l 
(2007) as a starting point to describe T(r), but this model was 
modified in two ways. First, we moved the sharp transition 
region (TR) down from a height z of 0.01 to 0.003 solar radii 
(Rq) to b etter match the cond i tions of semi-empirical mod- 
els (e.g.. iFontenla et al.l Il990t ICranmer & v an Ballegooiien 
120051; lAvrett & Loeseril2008l) . Thus, in the adopted model, at 
z = 0.01 Rq the temperature has risen to 0.48 MK, and it con- 
tinues to rise to a maximum value of 1.36 MK at z = 0.89/?q. 
We also increased the temperature slightly at distances greater 
than ^0.2 AU in order to better a gree with the mean o f the 
in situ T p and T e measurements of ICranmer et al.l (12009b . At 
r = 1 AU, T = 0.17 MK and it declines as T cx r' 6 . The 
one-fluid sound speed c s is defined as c 2 = ^k^T /;«h, where 
7 = 5/3 is the monatomic ratio of specific heats, k% is Boltz- 
mann's constant, and toh is the hydrogen atomic mass. 

Figure Q] shows the radial dependence of a selection of the 
background plasma properties defined above. It also shows 
the dimensionless plasma beta parameter, which is usually de- 
fined as the ratio of gas pressure to magnetic pressure, with 



1 gas 



mag 



V A 



(1) 



However, we will often use a simpler dimensionless parame- 
ter /? given by 



where f3 and /3q differ only by a factor of 1.2 when 7 = 5/3. 
The range of heights shown in Figure [TJ extends down into 



the solar chromosphere, but the wave models discussed below 
start at a lower boundary condition in the low corona; i.e., 
they specify the wave and turbulence properties only for z > 

O.O1.TV0. 

2.2. Linear Properties ofMHD Waves 

In this section we briefly summarize the dispersion proper- 
ties of linear MHD waves (i.e., phase and group speeds for the 
Alfven mode and the fast and slow magnetosonic modes) and 
the partitioning between fluctuations in kinetic, magnetic, and 
thermal energy. In Sections l2.3H2.4l we assume that all three 
types of MHD waves are present, and we vary their relative 
strengths arbitrarily in order to match the observations. 

The phase speed V p h = to/k is defined in terms of the fre- 
quency uj and the magnitude of the wavenumber k. In gen- 
eral, y p h is a function of the Alfven speed, the sound speed, 
and the angle 9 between the background field direction and the 
wavevector k. We follow the standard convention of defining 
a Cartesian coordinate system with the background magnetic 
field along the z axis and the k vector having components only 
in the x-z plane. Also, for now we express u and k in the frame 
comoving with the solar wind. For Alfven waves, 



v p 2 h 



vicos 2 e 



and for the magnetosonic modes, 



■(1±E) 



(3) 



(4) 



applies with the upper sign corresponding to the fast mode 
and the lower sign corresponding to the slow mode, and with 



£ = VL 



4/3 



(5) 



(1+/9) 2 

(see, e.g.. IWhang|[l997t iGoedbloed & Poedtsll2004h . In Sec- 
tion 12.31 we also need to know the component of an MHD 
wave's group velocity in the direction parallel to the back- 
ground magnetic field. We call this quantity V gz , and for the 
Alfven mode it is identically equal to Va no matter the value 
of 9. For the fast and slow modes, 



V ph cos6> 1 T 



a sin 9 
2S(1±S) 



(6) 



MHD waves excite oscillations in the plasma parameters. 
We denote the root-mean-square (rms) fluctuation amplitudes 
in velocity as v x ,v y ,v z , in magnetic field as B K1 B yi B z , and in 
density as dp. We ignore fluctuations in the electric field be- 
cause their contribution to the total energy density tends to be 
negligible when Va <C c. The kinetic, magnetic, and thermal 
energy densities associated with each type of fluctuation are 
given as 



PQVf 



Mi = 



8- 



e = /3 



(7) 



respectively, with i = x,y,z. For linear Alfven waves, the 
total energy density Ua is divided equally between trans- 
verse kinetic and magnetic fluctuations along the y axis, with 
U A = K y +M y and 

K L = M L= l_ 
U A ~ U A " 2 ■ 
For fast and slow mode waves, 



(8) 



Uf.s = K x +K z +M x +M z +e 



(9) 
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and we follow Whang (fl997h in expressing the partition frac- 
tions as follows, 



— — = /„ sin 2 9 + f t cos 2 9 , —±- = f„ cos 2 9 + /, sin 2 



fn 



c?A 



Mx 
Uks 



(V 2 h - C 2 )cos 2 ^ 
A 

e 

Uks 



y 2 (y p 2 h -c 2 ) cos 2 6 

(V p 2 h - C 2 )sin 2 < 



y^-yl 



(10) 

(11) 

(12) 
(13) 



where A = 4V 2 h - 



■2Vl- 



-2c 2 . The fast and slow velocity fluctu- 
ations (K x +K z ) always occupy exactly half of the total energy 
density, and the combination of magnetic and thermal fluctu- 
ations (M x +M z +Q) take up the other half. 

The energy partition fractions given above are fa miliar 
components of plasma physics an d MHD textbooks (e.g. JStixl 
11992k iGoedbloed & Poedtsll2004h . However, it is difficult to 
see intuitively how these fractions vary throughout the helio- 
sphere from Equations ([T0l>-([T3l alone. Thus, in Figure |2]we 
provide a schematic illustration of the energy partitioning for 
fast-mode waves. The three columns indicate the variation 
from low (/3 <C 1) to medium (/3 = 1) and high (fi 3> 1) beta 
plasmas. The three rows show the results for purely paral- 
lel propagation (9 = 0), an isotropic distribution of wavenum- 
ber vectors (see below), and purely perpendicular propagation 
(9 = tt/2). In general, all five terms on the right-hand side of 
of Equation (O are nonzero, but fractions less than ~ 1 % are 
not shown in Figure [2] This diagram can be transformed to 
show the properties of slow-mode waves by replacing (3 with 
1/(3 and interchanging the x and z subscripts with one another. 



2.3. Radial Transport Equations 

In order to determine how the total energy density of 
a given wave mode evolves with heliocentric distance, we 
solve equations of wave action conservation that contain mul- 
tiple sources of wave damping. There have been many 
discussions of energy conservation for both pure acous 
tic waves and incompressible Alfven waves (e.g., iDewat 

■ i ' i i i i i i i i > ~ 



1970; Isenberg & Hollwed 11982 IV elli 1993 UTu & Marschl 
1995MVerdini & VeilfeOol ISokolov et al.ll2009l) . but general 
derivations that can also be applied to fast and slow mode 
waves (for arbi trary 9) are less frequently seen. We utilize 
the results of Jacques] ( 1 19771) to write the damped wave action 
conservation equation as 



d_ 

dt 



tt) A dr 



{u + V gz . m )A U„ 



(14, 



where the subscript m can be replaced by A, F, or S for the 
relevant mode, Ao is the cross sectional area of the flux tube 
(i.e., Ao oc 1 /By), and Q m is the total dissipation rate for the 
mode in question. The dimensionless factor that takes account 
of the "stretching" effect of wavelengths in an accelerating 
reference frame is 
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FIG. 2. — Illustration of how fast-mode MHD waves divide their total fluc- 
tuation energy into kinetic, magnetic, and thermal energy in various regimes: 
wavevectors parallel to Bo (top row), an isotropic distribution of wavevectors 
(middle row), wavevectors perpendicular to Bo (bottom row); |3< 1 (left 
column), j3 = 1 (middle column), and /j > 1 (right column). Pl ot ted areas 
are proportional to the partition fractions given in Equations (10}-{T3). Ki- 
netic energy fractions are denoted by v x and v z , magnetic energy fractions are 
denoted by B x and B z , and the thermal energy fraction is denoted by 'th\ 

and the angle brackets denote a weighted average over all an- 
gles, 

[decern 

{f) = Jd9sm9 ' (16) 

where here we consider outward propagating waves with < 
9 < tt/2. The factor of cos 9 in Equation (fl3T l comes from the 
difference between the wave frequency in the Sun's reference 
frame (o>o) and the comoving-frame frequency (uj = ujq— k- Uo) 
that a ppears in the d efinition of the wave action; see Section 
III of I Jacquesl (1 1 977h . 1 Equation (PT"4l) implicitly assumes that 
ujq remains constant, but it does not require the specification 
of any given value of ljq. 

Our use of weighted averages over 9 is derived from the as- 
sumption that wave power is distributed isotropically in three- 
dimensional k space. In Appendix lAl we discuss the motiva- 
tions for assuming such an isotropic distribution of wavenum- 
ber vectors (specifically for the fast-mode waves). For Alfven 
waves, this assumption has no impact on solving Equation 
(fT4b . since the arguments of both angle-bra cketed quantitie s 
given above are independent of 9 (see also lHollweglll974T) . 
Thus, one obtains the same result for Alfven waves whether 
one assumes a single value of 9 or the isotropic distribution. 
For fast and slow mode MHD waves, some quantities de- 
pend strongly on 9 and others do not. For example, slow- 
mode waves in low-beta plasmas have values of the angle- 
dependent quantity uq + V gZim that are always nearly equal to 
uq + c s . Figure Q] shows the radial dependence of (ko + V^f) 
for the isotropic distribution of fast-mode waves. 

We solve Equation ( [Pfl ) for the energy densities of the three 
MHD modes (Ua, Up, Us), and we compute the dispersion and 

1 We note that the adopted form of Equations J 141 — 1 1 5t is only one out of 
several possible ways of placing and grouping the angle brackets. For the fast 
and slow modes, there is also potential ambiguity about whether one should 
use Lagrangian or Eulerian averages for uq in the transport equation. In future 
work we will explore the consequences of different methods of averaging. 



SOLAR WIND PROTON, ELECTRON, AND ION HEATING 



5 



energy partition properties of all three wave types as given in 
Section lZ2l At this stage, we neglect couplings between mul- 
tiple modes and other nonlinear effects. This is an approxima- 
tion that is likely to b reak down wherever the wave amplitudes 
become large ( e .g.. IChin & Wentzefl [T 972: Wentzel Il974t 
Goldstein! Tl 978: La combe & Mangenevl Il980t iPoedts et alJ 
1998t IVasquez & Hollwed 11999; De l Zanna et al.l 120011: 



Gogoberidze et al. 1 120071) . In Section g] we discuss the like- 



lihood of rapid coupling between the high-wavenumber tails 
of the Alfven and fast-mode power spectra. However, we con- 
tinue to assume that the total energy densities are given by the 
solution of the individual transport equations. 

To specify the dissipation rates Q m , we include both linear 
collisional effects (e.g., viscosity, thermal conductivity, and 
electrical resistivity) for all three modes and nonlinear turbu- 
lent damping for the Alfven and fast mode. Thus, we use 

Qa = Qa + 2iaU k , Q F = Q F + 2~f F U F , Q s = 2 ls U s . (17) 

We give the amplitude damping rates j m , which include an 
approximation for the transition from strongly collisional to 
collisionless regimes, in Appendix |B] The turbulent damp- 
ing rates (5a and Q F are described in more detail below. In 
general, these rates depend on the parallel and perpendicular 
components of the wavenumber (k\\,k±)- For the purposes of 
evaluating these rates in the global wave transport equations, 
we assumed that k± = 1 / X± for all three modes, where A^ is 
the turbulent correlation length described below. For the fast 
and slow modes, our assumption of an isotropic distribution of 
wavenumbers is consistent with also assuming = k±. For 
the Alfven mode, we found that 7a never depended on the 
assumed value of k\\ at all, but for completeness we used the 
critical balance condition (introduced in Appendix© to spec- 
ify k l{ . 

We adopt phenomenological forms for the turbulent dis- 
sipation rates that are equivalent to the total energy fluxes 
that cascade from large to small scales. Thus, (5a and 
(5f are constrained only by the properties of fluctuations at 
the largest scales, and they do not specify the exact kinetic 
means of dissipation once the energy reaches the smallest 
scales (but see, however, Section |5). Dimensionally, these 
are similar to the rate of casc ading energy flux derived by 
Ivon Karman & Howarthl d 1 9381) for isotropic hydrodynamic 
turbulence. For the nonlinear dissipation of Alfvenic fluctua- 
tions, we use 



Qa = PO&A^tui-b 



zlz++zlz. 

4Al 



(18) 



(see also iHossain et alJ 119951: IZhou & Matfhaeusl [l990a; 
Matthaeus et"afl 119991: iDmitruk etakl 1200 ll 12002 



Breech et al . 2008). For the fast-mode waves, we use 



. (v, 2 + v?) 2 



(19) 



where the quantity (v^ + vf) collects togethe r the total kinetic 
energy in fast-mo de velocity fluctuations (Chandran 120051; 
Suzuki et al. 2007). Many of the terms introduced in Equa- 
tions (|L8T>— (fl~9b are defined throughout the remainder of this 
subsection. 

Equation (fT~8~b depends on the magnitudes of the [Elsasser 
(1950) variables, Z± = v y ±B y /(4irp ) l/z , which specify the 
power in outward (Z_) and inward (Z + ) propagating Alfvenic 
fluctuations. Alfvenic turbulent heating occurs only when 
there is energy in both modes. In practice we compute an ef- 
fective reflection coefficient TZ = |Z+|/|Z_| whose magnitude 



is always less than unity, and thus we express the Elsasser 
variables in terms of the Alfvenic energy density as 



Po(l+^ 2 ) ' 



TZZ- 



(20) 



An accurate solution for Z± requires the integration 
of non-WKB equations o f Alfven wave reflec tion (e.g., 
iHeinemann & Olbertl [19801 iVerdini & Veflil 120071) . How- 
ever, our assumption that the total power Ua varies 
in accord with straightforward wave action conserva- 
tion has been shown to be reasonable, even in envi- 
ronme nts where TZ is not small such as the chromo- 
sphere (van Ballegooiien et al.l201 ll) and interplaneta ry space 
dZank et alJll996tlCranmer & van Ballegooiienll2005l) . 

We estimate the reflection coefficient TZ using a 
modification of th e low -frequency approximation of 
IChandran & Hollwegl (120091) . Specifically, we examine the 
magnitudes of terms in the transport equation for the inward 
Elsasser variable, 

dZ + dZ+ f Z+ Z_ \ Z+Z_ 

_^ + (M0 _y A )_i = (M0 + y A )^ + _j-_l_ (21) 



where 



V A 



oVA/or opo/or 

IChandran & Hollwegl (120091) neglected both terms on the left- 
hand side of Equation (fJTJ as well as the term containing Ho, 
and thus were able to solve for Z + straightforwardly. However, 
in cases of strong reflection, the term containing H D may have 
a magnitude comparable to the other dominant terms. Thus, 
we keep all three terms on the right-hand side and solve for 



where 



n 



h 



2h/\H A \ 
J l + (h/\H D \) 

A_l(mq + V a ) 
2Z_ 



(23) 



(24) 



Equation (22) of Chandr an & Hollwegl (120091) is recovered in 
the limit of h <C I Hp I, with TZ » 2 h /\Ha | . In the case of purely 
linear reflection. ICranmerl (120 1 Ol) found that the most accurate 
local estimates for TZ were obtained when Ha was replaced 
with the positive-definite quantity 



#A = VAfref = (r + R G ) 1 



(25) 



We used H A instead of H A in Equations d23l-(l2"4"ii to compute 
TZ. 

The definitions of the turbulent dissipation rates contain the 
perpendicular length scale Aj_, which is an effective trans- 
verse correlation length of the turbulence for the largest "outer 
scale" eddies. For simplicity we use the same correlation 
length for both the Alfvenic and fast- mode fluctuations, but 
this may not be universally valid (e.g., Suzuki et al. 2 0071) . In 
previous papers we assumed that Aj^ scales with the trans- 

-1/2 

verse width of the magnetic flux tube; i.e., that X± on B Q 
(lHollwegl ll986l). Here we describe the evolution of the trans- 
verse correlation length A^ with the following transport equa- 
tion, 



d\ ± Aj_ dA j3 A (ZlZ^Z\Z. 



dr 2A dr u + V A \ Z^+Zl 



(26) 
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where /3a is a dimensionl ess constant that is often assumed to 
be equal to a\/2 (e.g.. iHossain et aDl!995h . The first term 
on the right-hand side of Equation (l26l i drives the correla- 
tion length to expand linearly with the perpendicular flux-tube 
cross section (lHollweglll986l) . The second term takes account 
of the nonlinear coupling between the fluctuations and the 
background plasma properties. It is given in a form suggested 
initially by Mattha eus et al.l d 19941) and later generalize d to 
nonzero cross-helicity turbulence by iBreech et al.l (12008b and 
others. Our transport equation attempts to bridge together the 
effects of the two terms. In the lower solar atmosphere (be- 
tween the photosphere and the chosen lower boundary of 0.01 
Rq for the wave transport models) we assumed that the first 

1/2 



term in Equation ( 1261 1 is dominant, and thus Ax oc A, 

The turbulent dissipation rates also depend on dimen- 
sionless Kolmogorov-type constants a& and a? that are of- 
ten assumed to hav e v alues of order unity . For example, 
IHossain et all d 19951) and lBreech elaT] ( 120091) found that a A « 
0.5 gives rise to dissipation rates that agree well with both 
numerical simulations and heliospheric observations. In our 
case, we used this value as a starting point, but we also var- 
ied qa as a free parameter in order to produce the best match 
to the well-constrained Alfvenic fluctuations. On the other 
hand, the properties of heliospheric fast-mode turbulence are 
not known nearly as well as the Alfven- wave turbulence. We 
thus relied on the independent wave-k inetic simulations of 
Pongkitiwa nichakul & Chandranl (1201 2l) to fix a-p at a value 
of 2.3. 

The Alfvenic cascade rate contains an efficiency factor £ tur b 
that attempts to account for regions where the turbulent cas- 
cade may not have time to d evelop before the fluct uations are 
carried away by the wind. Cranmeretal. (2007) estimated 
this efficiency factor to scale as 



1 



turb 



1 + (/eddy/ ftef) 



(27) 



where the two timescales above are f e ddy, a nonlinear eddy cas- 
cade time, and f re f , a timescale for large-scale Alfven wave re- 
flectio n (see also Dmitruk & Matthaeus 2003; Oughto net al.l 
2006). The reflection time is often defined as f re f = 1 /|V • \a\, 
but we solved Equation ( f25l ) for f re f in order to remain consis- 
tent with the adopted model for 1Z. The eddy cascade time is 
given by 

feddy = u+MaK ' (28) 

where the Alfven Mach number M A = wo/Va and the nu- 
merical factor of 3ir comes from the normalization of an as- 
sumed shape of the turbulence spe ctrum (see Appendix C of 
ICranmer & van Ballegooiien2005). The two limiting cases of 
£turb •C 1 and £ t urb ~ 1 are roughly equivalent to the "weak" 
and "strong" cascade phenomenologies discussed in Section 
[3] but they are not precisely the same. 

2.4. Representative Solutions 

We solved the transport equations given in Section |2~31 by 
numerically integrating upwards from a specified set of lower 
boundary conditions at z = 0.01 Rq and assuming time-steady 
conditions (i.e., dU,„/dt = 0). We used a logarithmic grid of 
500 radial zones in z that expands out to a maximum distance 
of 86O/?0 « 4 AU. The transport equations were solved with 
straightforward first-order Euler steps. The values of the El- 
sasser variables Z± in each zone were determined by iteration, 



TABLE 1 

Standard Model Parameters for 
Coronal Hole MHD Wave 
Transport 



Parameter 


Value 


&A 


0.60 


h 


0.31 


dip 


2.3 


(f/ A /po) 1/2 (atz = 0.01R Q ) 


29.0 km s _1 


(t/p/po) 1 / 2 (atz = 0.01R Q ) 


24.3 km s -1 


(U s /p ) [/2 (atz = 0.01R©) 


9.17 km s- 1 


Ax (at photosphere) 


120 km 



since Equations ( f20b and (l23l do not give a simple closed- 
form solution for Z+ and Z_ by themselves. 

There are a number of free parameters in this model whose 
values were not easily obtained from either theoretical cal- 
culations or observations. In addition to the lower boundary 
conditions on the wave energy densities Ua, Up, and t/s, there 
is also the lower boundary condition on the correlation length 
Ax and the values of the two von Karman constants a\ and 
$a- Initially, we varied these six parameters randomly in or- 
der to build up a large Monte Carlo ensemble of trial solu- 
tions. For each model, we synthesized the radial variation of 
observable plasma fluctuations such as the root mean squared 
(rms) parallel and perpendicular fluctuation speeds, 



vii = v. 



1/2 



(29) 



the Elsasser variables Z±, and the rms fractional density fluc- 
tuation amplitude Sp/po. The velocity amplitudes vii and vx 
contain contributions from all three MHD wave types. In 
nearly all models produced here, vii is dominated by the fast 
mode and vj_ is dominated by the Alfven mode. Observations 
of these quantities are discussed below. 

There was no single set of parameter values that gave rise 
to perfect agreement between all of the synthesized and ob- 
served fluctuation quantities. This is not surprising, since 
the models are certainly incomplete and there are signifi- 
cant uncertainties in the observations and their interpretation. 
Also, even though we aimed to restrict ourselves to measure- 
ments made in "quiet" high-latitude fast wind streams, some- 
times only low-latitude data were available. Thus, in Table 
Q] we give a set of optimized parameters that were chosen 
because they produce adequate agreement with the full set 
of observed quantities. There were other combinations of 
the six parameters that gave better agreement on any single 
observation, but in most of these cases the agreement be- 
came worse for other observations. Although the ratio of 
the two von Karman constants &a/$a was allowed to vary 
freely, the optimal value was non etheless found to be close 
to the commonly used value of 2 (IHossain et al]ll995l) . The 
best photospheric value of Ax ~ 120 km is inte rmediate be- 
tween the values of 75 km (ICranmer et a l. 2007) and 300 km 
(ICranmer & van Ballegooijen 2005) found from earlier mod- 
els. 

Figure [3] shows the comparison between synthesized and 
observed fluctuation quantities for the model parameters 
given in Table Q] The observational constraints on vj_ at 
z < O.I^q are a combination of the off-limb nonth ermal 
emission line widths g iven by IBanerjee et alJ (119981) and 
ILandi & Cranmerl (|2009|). The ob s ervat ions shown between 
0.3 and 1 Rq are from Esser et al.l (1 19991) . At larger heights, 
vx becomes approximately equal to Z_/2, so we truncate the 



SOLAR WIND PROTON, ELECTRON, AND ION HEATING 



7 




0.01 0.1 1 10 100 1000 

(r/R ) - 1 



1 

_(b) 


i i 


i 


1 


I 






S y 
■ N / 
\\/ 

^ - v 

"\ N 
\ 


\ \ 

V \ \ 




i 


(cm) 

i i 


i 


A> \ 

\\ \ 

% X 
\\\ 

\\\ 

\\ 

i 


I 



0.01 0.1 1 10 100 1000 

(r/R ) - 1 

FIG. 3. — (a) Model values for Vj_ (black solid curve), Z- (black dot- 
dashed curve), and Z+ (black dotted curve) compared with measurements 
(gray boxes). Model values for vn (red dashed curve) compared with mea- 
surements (light red circles). Velocities are plotted in units of km s ; see 
Equation (29). Total density amplitude Sp/po (blue solid curve) is shown 
with its components from fast-mode (blue dot-dashed curve) and slow-mode 
(blue dotted curve) waves, and compared with observations (light blue re- 
gions and rectangles), (b) Modeled total heating rate 2toi/po (red dashed 
curve) compared with empirical constr aints (light red regions) and the total 
heating rate from Cranmer et al. (2007 ) (green dot-dashed curve), all in erg 

s g -1 . Standard model value for Ax (solid black curve) compared with ear- 

-1/2 

lier assumption Ax oc B (black dotted curve) and with in situ estimates 
(gray region), shown in units of cm. See text for data sources. 

vj_ curve in favor of showing the radial dependence of Z+ and 
Z_ more clearly. The latter are compared di rectly with high- 
speed wind data from Helios and Ulysses dBavassano et alj 
2000). Observations of longitudinal velocity fluctuations are 
more difficult to find, and we show only t he o n-disk nonther- 
mal line width velocities of IChae et all (1 1998b as a way to 
compare with the modeled values of vii . 

Figure[2a) shows how the modeled density fluctuation am- 
plitude 5p/po is dominated by slow-mode waves in the low 
corona (z < OARq) and by fast-mode waves in the extended 
corona and solar wind (z > IRq)- The low-corona observa- 
tions are drawn as an approxi mate boundary reg ion around 
the polar plume data given by lOfman et al.l (11999b . The in- 
termediate data point at z = 4-Rq is an e mpirical value of 
Sp/po estimated from radio sounding data (IColes & Ha rmon 



Il MlSpanderl2002tlHarmon & Co les 2005; IChandran et al.l 
2009), but it is still unclear what fraction of the measured den- 
sity fluctuations are due to anything even close to ideal MHD 
waves. At larger distances, we sh ow approximate ran ges of 
density fluctuations as reported by Marsch & Tu (1990) (blue 
rectangl es at z < 200/^1 iTu & Marschl (11994 (open rectan- 
gle), and llssautier et a l. ( 1998) (blue rectangle at z > 300R Q ). 
FigurefJJb) compares the result of solving Equation (l26l for 

-l 11 

Aj_ with the simpler approximation of Ax oc B Q 1 . The plot 
also shows fa st-wind estimates o f Aj_ between 1 .4 and 5 AU 
from Ulysses dBreech et al.ll2008l) . Figure |3jb) also compares 
the total heating rate Q (ot = Qa + Qf + Qs with observational 
constraints and with the modeled coronal heating rate from 
ICranmer et al.l (120071) . The shaded area between 0.2 and 5 Rq 
is an envelope surrounding a collec t ion of empirical and the- 
oretica l heat ing curves f rom Wang (11994b . lHansteen & Leerl 
d 19951) . and lAllen et al l (119981) . These rates illustrate what 
is needed to produce the observed coronal heating and so- 
lar wind acceleration. The area shown at larger distances 
(z > 6QRq) is a representation of the range of to tal (proton and 
electro n) empirical heating rates estimated by Cran mer et al.l 
(2009). Note that the turbulent heating rates Qa and Qp dom- 
inate the total heating rate, with approximately 70% of the 
total coming from Qa and 20% from Qp. Less than 10% of 
2 to t comes from the linear damping terms. 

There are additional measurement techniques that may be 
used to further constrain the model parameters, and in future 
work we will incorporat e as m any of these as possible. For 
example, Holl weg et al.l (2010) argued that radio measure- 
ments of Faraday rotation fluctuations may put unique em- 
pirical constra i nts on the value of Ax in the corona. Also, 
Sahra oui et al.l (('2010) used multi-spacecraft data to tease out 
new details of the wavenumber anisotropy of MHD fluctu- 
ations, which may lead to better limits on, e.g., Vm/v_l in 
the heliosphere. Unfortunately, the vast majority of these 
measurements have been made for the slow solar wind and 
not the much less structured fast win d associated w i th po- 
lar coronal holes. Nearer to the Sun, Kita gawa et al.l (1201 Ol) 
used the dispersive and energy partition properties of thin- 
tube MHD waves to diagnose the presence and strengths of 
various modes in active regions. These techniques may be 
useful in open-field regions as well. 

Although we did not include any explici t mu lti-mode cou- 
pling in the transport equations of Section 12.31 there is some 
feedback between the modes. For example, the correlation 
length Ax is used in both the Alfvenic and fast-mode turbulent 
heating expressions, and it is also used to set the wavenumbers 
and k± in the linear dissipation rates j m . Thus, the choice 
of the lower boundary condition on Ax can have a significant 
impact on the radial evolution of all three wave types. Figure 
2] illustrates this by varying the photospheric value of Aj_ be- 
tween 30 and 300 km and using the other standard parameters 
from Table [TJ The integrated energy densities are plotted in 
velocity units as (U m / po) 1 ^ 2 ■ The power in the Alfven waves 
changes by only a small amount because the damping is never 
a strong contributor to the Ua transport equation. However, 
damping is a major effect for the fast and slow modes, and 
thus small changes in the damping rate's normalization can 
have large relative impacts on the resulting energy densities. 

Figure [4] shows that, no matter the choice of normalization 
for Ax, it seems unlikely for the slow-mode waves to have 
significantly large amplitudes anywhere but in the lowest few 
tenths of a solar radius. This appears to be consistent with 
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FIG. 4. — Radial dependence of MHD wave energy densities per unit mass, 
for a range of photospheric boundary conditions on on Ax- From bottom to 
top in each set of curves, the values are: 30 km (black), 60 km (dark blue), 100 
km (cyan), 120 km (green), 200 km (orange), and 300 km (red). Different line 
styles denote Alfven waves (solid curves), fast-mode waves (dashed curves), 
and slow-mode waves (dotted curves surrounded by gray background). 

models of slow-mode shock fo rmation and dissipation in po- 
larplumes ( ICuntz &~S uess 2001). Therefore, in the remainder 
of this paper, our models of turbulence in the fast solar wind 
ignore the slow-mode waves altogether. We also note that Fig- 
ure|4]suggests that the actual fast-mode wave properties in the 
high-speed solar wind may be more highly variable than the 
Alfven wave properties. Our use of a "standard" model for 
the fast-mode waves (using the parameters given in Table Q3 
is thus presented as an example case and not a definitive pre- 
diction. 

3. MHD TURBULENT CASCADE 

In this section we begin constructing a model of the 
wavenumber distribution of Alfven and fast-mode fluctuation 
power at each radial distance. We make use of a general as- 
sumption of "scale separation;" i.e., we presume that the tur- 
bulence becomes fully developed on timescales short com- 
pared to the bulk solar wind outflow and the large-scale ex- 
pansion of open flux tubes. This allows us to model the tur- 
bulence as spatially homogeneous in a small volume element 
with constant background plasma properties. This seems to 
be the general assumption made by the majority of MHD 
simulations of turbulence in the solar wind. 2 Whether or not 
this approximation is valid, it is useful to begin studying the 
wavenumber dependence of the cascade in this manner. 

3.1. Wavenumber Advection-Diffusion Equations 

We model the MHD fluctuations as time-steady Fourier dis- 
tributions of wave power in three-dimensional wavenumber 



2 See, however, "expanding box" type simulations (GrapDin & Velhl l 199(3 : 
Liewer et al. 200 1 ) that attempt to include some aspects of the large-scale 
radial evolution of the plasma parcel undergoing a turbulent cascade, and 
collisionless kinetic models that include expa nsion effects together with local 
diffusion in velocity space (Isenberg & Vasquez 2009, 2011). 



space. Although additional information about the physics of 
turbulence can be found in more complex statistical measures 
of the system (e.g., higher-order structure functions), we limit 
ourselves to describing the power spectrum because that is 
the basic quantity needed to compute the quasilinear particle 
heating rates. 

Because of the simplified flux-tube geometry discussed in 
Section im we assume the background magnetic field is par- 
allel to the bulk flow velocity, and thus the syste m has only 
one p referred spatial direction (see, however, iNarita et all 
2010). The random turbulent motions create a statistical 
equivalence between the x and y directions transverse to the 
background field, so that we can describe the power spectra 
as two-dimensional functions of kn and k± only. By conven- 
tion, we define the full three-dimensional power spectrum E„, 
in effective velocity-squared units; i.e., when integrated over 
the full volume of wavenumber space, the spectrum gives the 
fluctuation energy density per unit mass, or 



— = /of 3 k£ m (k). 

Po 



(30) 



In Appendix [A] we review some of the basic physical pro- 
cesses that determine the shape of the spectrum for Alfvenic 
(m = A) and fast-mode (m = F) fluctuations. 

We describe the driven turbulent cascade as a combina- 
tion of advection and diffusion in wavenumber space. At 
first, it may appear that a smooth and continuous descrip- 
tion of the spectral "spreading" of a cascade ignores too 
much of the inherentl y stochastic and non local nature of 
turbulence. However, Chandrase kharl (1 19431) showed that 
such a model can be made to capture the essential statis- 
tics of a large ensemble of random-walk-like (i.e., Brow- 
nian) processes. Specific models of turbulent wavenum- 
ber transport using diffusion or advecti on equati ons in- 
clude thos e of IPaoT (fl%5l). iLeithl ([19671). iTu et al] ([19811) . 
Tul d 19881). IZhou & Matthaeusl dl990bl). iMiller et al] dl99fj|) 
Stawicki et alJ (1200 ll). IChandranl (l2008bl), Matthae usetatl 
(120091) . Uiang et all (120091) . and lGaltier & Buchlinl (l2010h . For 
the cascade of Alf venic fluctuations, we general l y follo w the 
approach taken by Cran mer & van Ballegooiien (2003). The 
general forms of these equations are given as 



dE A 
dt 



1 d 



d 



k± dk± 



(k 2 ± E A 



-/'J 



+ a 



dk n 



D 



All 



dEp 
dk\\ 

dEp 
~dk 



+ Sa-2 7a £ a + C af (31) 



and the terms on the right-hand sides of Equations (|3H-(|32l 
are defined throughout the remainder of this subsection. The 
mode coupling term Caf is described further in Section |4] and 
the dissipation rates 7a and 7f are described in Section [5] 

The perpendicular Alfvenic cascade is described by the first 
term on the right-hand side of Equation d3lV and we as- 
sume a n arbitrary linear combination of a dvection and dif- 
fusion. ICranmer & van B allegooijenl (120031) found that many 
key properties of the turbulence do not depend on whether the 
cascade is modeled as advection, diffusion, or both, so we re- 
tain all terms for maximum generality. For both the parallel 
Alfvenic spectral transport and the isotropic fast-mode trans- 
port, a more standard diffusion coefficient is assumed. The 
dimensionless multipliers to the E A diffusion coefficients are 
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denoted a± and au, to correspond roughly to 5a in Equation 
( fT~8b , and the dimensionless multiplier for the wavenumber 
advection coefficient is denoted 

For the Alfvenic cascade, the overall behavior of wavenum- 
ber transport in the perpendicular and parallel directions is 
specified by the diffusion-like coefficients 



TA 



D A || = 



V A 



(33) 



where ta is the cascade timescale defined below, and v± is the 
A:^ -dependent velocity response of the waves. Note that D A || 
is independent of s o it can be pulled out of the deriva- 
tive in Equation ( Bil l. ICranmer & van Ballegooiie"nl (12003b 
showed that the above form for the diffusion c oefficients tends 
to r eproduce the IGoldreich & Sridharl (fl995 ) critical balance, 
and iMatthaeus et alj d2009f) derived similar functional forms 
for the coefficients. When specifying the properties of the 
wavenumber cascade, we apply the scalings for "balanced" 
turbulence (i.e., zero cross helicity, or Z+ = Z_), which is more 
straightforward to implement but is formally inconsistent with 
the large-scale transport model of Section |2] 

For ideal MHD Alfvenic fluctuations, v\ is equal to b\, the 
latter representing the transverse magnetic variance spectrum 
divided by Airps to convert it to units of velocity squared. Fol- 
lowing the usual convention, the power spectrum E A tracks 
the magnetic fluctuations, so the reduced spectra are defined 
formally as 



b\ = k\ 



(34) 



The dimensionless factor 4> describes the departure from ideal 
MHD energy equipartition. For small values of k±, we as- 
sume <j> m 1. However, as k± increases into the regime of 
kinetic Alfven waves (K AWs), <fi can become much larger 
than 1 . iHollwejd ( 1 1999b described how the main difference 
between v_l and b± in the KAW regime comes from an en- 
hanced response of the electron velocity distribution to the 
electric and magnetic fluctuations. For simplicity, we use an 
approximate analytic expression 



k 2 V 2 
II A 



\+k\plm e /{pm p ) ' 



(35) 



where p p = w p /Q, p is the proton thermal gyroradius, with the 
proton most-probable speed given by w p = {2k^T p /m p ) l l 2 and 
the proton cyclotron frequency by Q, p = eBj m p c. Our term tj> 
is equivalent to a 2 as defined bv lHowes etail d2008l) . 

Inspired by Equation (IA6t . we define the Alfvenic spectral 
transport timescale as 



ta = 



kjyj_ 



(36) 



where we chose to replace the general critical balance param- 
eter x by its value at the outer-scale parallel wavenumber & || ■ 
Thus, 

w o k \\V A 
Xo = t— ~ tt- , (37) 



k±v_ 



k\b\ 



and xo i s me appropriate critical balance parameter to 
use when solving for the properties of the dominant low- 
frequency cascade. From Equation ( |35l l we see that KAW 
outer-scale frequency ujq I//2 £ O || V\, so that a factor of </>'/ 2 
cancels out of both the numerator and denominator to give 



the final approximate expression above. The wavenumber &o|| 
specifies the spatial scale along the field at which energy is 
injected in the source term 5a (see below). Because fe || is 
assumed to be constant (at a given heliocentric distance r), 
the parameters xo and ta are both functions of k± and not 
k\\. The ab ove form for Equati o n d36l was motivate d by th e 
ana lysis of IZhou & Matthaeusl d!990bl) . IChandranl d2008bl) . 
and lHowes et al.l d2012l) . who described how the cascade and 
wavenumber anisotropy change when the system transitions 
from weak (xo ^> 1) to strong (xo *C 1) turbulence. 

As mentioned above, our expressions for t a , D a ±, and 
Da|| assume zero cross helicity (i.e., 71= 1). There is still 
no agreement about how to generalize these terms when in- 
efficient wave r eflect ion gives rise to nonzero cross helicity. 
iLifhwick et al.1 d2007l) found that the cascade timescales for 
outward and inward wave modes are different from one an- 
other when 1Z ^ 1, but their parallel spatial scales are the 
same. However, Beresnvak & Lazarian (2008, 2009) found 
that for the outward mode s hould be larger than for 
the inward mode, and thus the IGoldreich & Sridharl (1995) 
critica l balance must b e modified (see Equation (l45l l be- 
low). IChand ran (2008b) outlined a method for setting up the 
advection-diffusion equations in the case of 1Z 1, but we de- 
fer a full implementation of that approach to future work. 

Putting aside the issue of imbalanced turbulence, the dom- 
inant perpendicular nature of the Alfvenic cascade allows us 
to define a reduced transport equation that follows the evolu- 
tion of the spectrum as a function of k± only. If we ignore the 
mode coupling term Caf for now, we can multiply Equation 
( f3TT > by k\ and integrate over k^ to obtain 



db\ 

dt 



d 
dk~ 



db\ 
at i 



-H±bj 



+ S A -2j A b A 



(38) 

This is essentially the s ame as Equation (11) of 
Cranmer & van Ballegooiien (2003). The reduced source 
term S A and dissipation rate 7a are defined similarly to the 
corresponding terms in Equation (1311 , but they are weighted 
toward the low-£|| regions of wavenu mber spac e that are 
"filled" by the cascade. In Appendices IC.H4C.3I we derive 
analytic solutions for the time-steady Alfven-wave power 
spectrum in various limiting cases. 

The cascade of fast-mode waves, described by Equation 
(|32| >. appears to be conceptually simpler than the strongly 
anisotropic Alfven-wave cascade. The diffusion coefficient is 
given by Dp = k 2 /rp, wh ere r p is related to the IK-like cascade 
time given by Equation dA2t with p = 1 . There is increasing 
evidence (e.g.. iMarkovskii et al.ll2.Q10l) that a fast-mode cas- 
cade is more rapid in the directions perpendicular to the field 
than along the field. However, the cascade does appear to pro- 
ceed outward "radially" in the direction of increasing k. Thus, 
it makes the most sense to use an isotropic diffusion formal- 
ism as in Equation d321 l. but scale the magnitude of the diffu- 
si on timescale with 0. Following the weak turbulence model 
of IChandranl (120051) . we adopt 



7"F 



Va 



kvj sin 9 



(39) 



fcVsintf 4mt 6 £ , F sin( 



which implies that 

D F = = 

Va Va 

IChandranl (120051) showed that the sin 9 dependence in the de- 
nominator of Tp is consistent with an isotropic energy flux for 



(40) 
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the cascade, but it does not guarantee an isotropic wavenum- 
ber spectrum Ep(k). More information about how we cho se to 
implement the fast-mode cascade is given in Appendices IC. 41 
andlC31 

In order to fully describe the cascade in the advection- 
diffusion equations, four dimensionless spectral trans- 
port constants (ctj_, M-L> a F) need to be specified. 
Matth aeus et alj d2009l) summarized the results of many MHD 
turbulence models and found that a± often takes on values 
between 0.2 and 0.5, and au « 0.43 a_i_ seems to be a usefu l 
parameterization (see Equa tion 13 of Matth aeus et al.l|2009l) . 
IZhou & Matthaeus! d 1990b!) and lMatthaeus et alj (120091) made 
a case for a classical form of the diffusion operator th at im- 
plies = 2a±. Alternately, van Ballegooiien (1986) found 
that a cascade of random-walk-like displa cements of m a gnetic 
flu x tubes is describe d well by fi± = a±. iHowes et al.l (2008) 
and lChandranl d2008bl) used a straightforward advection equa- 
tion to model an Alfvenic cascade, whic h sets a± = and as- 
sumes (j,± y 0. For this type of model, IHowes et alj (120081) 
derived /i± w 0.2. 

In our models, we are constrained by the values of the cas- 
cade constants a& and dtp used in the global transport equa- 
tions of Section [2] We related these constants to the ones 
defined above by integrating the cascade advection-diffusion 
terms over wavenumber to find dll m /dt. By demanding this 
quantity be equal to the heating rate Q m , we obtained 



2a± 



■ + A 1 j 



3V6tt 



(41) 



which assumes that the perpendicular cascade is dominant and 
that 7?. ss 1, and 

32 



OtV = ■=- «F 
V7T 



(42) 



We keep the ratio s = fi±/ a± as a free parameter and we 
explore the ramifications of varying it b elow. Note, how- 
ever, th at if we used s = 2 (as assumed by IZhou & Matthaeus! 
1990b), then Equation (gT]l gives a± 0.73 and n± « 1.47. 
These are roughly consist ent with the constan t s give n by 
IZhou & Matthaeus! d!990bl) and iMatthaeus etail d2009l) . To 
complete the system of cascade constants, we adopt the 
IMatthaeus et~aTTd2009l) choi ce for ot \\ = 0.43qj_, bu t we com- 
pute this quantity using the IMatthaeus et al. (2009) assump- 
tion of s = 2. 

The source terms, S\ in Equation (|3"TT i and Sp in Equation 
d32l . describe the outer-scale injection of fluctuation energy. 
The global energy balance of the waves is already described 
by the radial transport model of Section |2] Thus, we spec- 
ify the magnitudes of Sa and Sp by demanding that the time- 
steady total energy densities Ua and Up be maintained at their 
known values at a given distance r. From a physical stand- 
point, however, it is unclear whether the passive propagation 
of waves dominates the source terms, or whether there is sig- 
nificant local "stirring" that converts large-scale dynamical 
motions (e.g., velocity shears in evolved corotating streams) 
into new fluctuations. 

We adopt specific functional forms for S\(ku , kx) and Sp{k) 
that are described in detail in Appendix [C] Generally, the 
source terms are nonzero only at the lowest wavenumbers, at 
which the fluctuations are driven. For the Alfven waves, we 
continue to use the assumption from Section l2~3l that the per- 
pendicular driving scale is set by the turbulence correlation 
length; i.e., ko± = 1/Aj_. For the fast-mode fluctuations, we 
assume their outer-scale wavenumber magnitude £qf is also 



equal to ko±, since the largest-scale transverse stirring mo- 
tions are likely to be common to both Alfvenic and fast-mode 
waves. There are several ways that one could imagine defin- 
ing the parallel outer-scale Alfven wavenumber fc ii : 

1 . Monochromatic Alfven waves that propagate up from 
the corona retain a constant frequency ojq in the Sun's 
inertial frame. However, because the phase speed varies 
with distance, the corresponding wavelength undergoes 
"stretching" commensurate with the dispersion relation 



ko\\ 



u + Va 



(43) 



2. The fluctuations propagating up from the Sun 
may already be fully turbulent (see, e.g., 
Ivan Ballegooiien et al.l 1201 ll) . Thus, the outer-scale 
parallel wavenumber may be coupled continuously 
to the perpendicular waven umber via critical balance 
dGoldreich & Sridhail995l) . with 



kpx 
V A 



(44) 




3. In flux tubes with nonze r o cross heli city (i.e., 1Z < 
1), iBeresnvak & Lazarianl (120081 120091) found that the 
inward waves should obey the iGoldreich & Sri dhar 
d 19951) critical balance, but the outward waves (which 
are generally what we intend to model) obey a modi- 
fied version of critical balance, which we approximate 



(45) 



4. In some cases we assume that the dimensionless ra- 
tio £o||Ao_i_ remains fixed at a constant specified value. 
Many studies of MHD turbulence assume isotropic 
forcing at the outer scale, which is consistent with the 
fixed ratio k 1| /ko± = 1. The lack of a physical justifica- 
tion for this approximation is offset by its simplicity. 

Figure [5] illustrates the ratio fco||Ao_L for several of the above 
methods of setting the parallel outer scale. For example, it 
shows the result of evaluating Equation (1431 for a range of 
wave periods P = 27r/wo between 1 and 100 minutes. Constant 
assumed values of k { )\\/ko± would correspond to horizontal 
lines in Figure [5] 

3.2. Solutions in the Absence of Coupling 

Here we present some example results for the power spectra 
£a(£||,£±) and Ep{k^,k±). These spectra are computed from 
Equations d31t . d32l , and (1381 in the limiting cases of time 
independence and no mode coupling {Cap = 0). The Alfvenic 
spectrum was first computed in its redu ced f or m usin g the so- 
lutions for b±(k±) given in Appendices lC. 1 l and lC.2l and then 
it was expanded into full wavenumber space by using the re- 
sults of Appendix IC. 3 1 The shape of the fast-mode spectrum 
was determined from the analytic solutions given in Appen- 
dices |C3 and |C5] 

To illustrate the wavenumber dependence of the power 
spectra, we chose a single coronal height z = IORq at which 
/3 « 0.04. We typically plot the wavenumbers in terms of di- 
mensionless quantities £||Va/O p and k±p p . Dissipative wave- 
particle interactions tend to become important when these 
quantities reach order-unity values, and ideal MHD condi- 
tions apply when these quantities are small. Typically, the 
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FIG. 5. — Radial dependence of the modeled ratio of outer-scale 
wavenumbers &o||/Aox computed under various assumptions: constant 
inertial-frame frequencies (red solid curves, labeled by wave period), ideal 
IGoldrei ch & Sridhar ( 1995) crit ical balanc e (dotted black curve), and mod- 
ified Beresnyak & Lazarian (2008, 2009) critical balance (dashed black 
curve). 
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FIG. 6. — Reduced Alfvenic fluctuation spectra for magnetic field and ve- 
locity fluctuations at z = lOi?0, plotted as a function of k±p p . Undamped 
spectra for b± (red dashed curve) and v± (red dot-dashed curve) are com- 
pared with damped spectra for b± (black solid curve) and Vx (black dotted 
curve). The dimensionless KAW dissipation rate j/w used to compute the 
damped spectra is also shown (green solid curve), as is the location of the 
perpendicular outer scale &oxPp (blue dotted line). 

driving scale for Alfvenic turbulence occurs at koxPp ~ 10~ 6 
to 10~ 4 , with the larger values generally occurring at larger 
heliocentric distances. 

In Figure|6]we show the time-steady k± dependence for the 
Alfvenic b± and vx fluctuations, both with and without KAW 
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FIG. 7. — Reduced Alfvenic magnetic spectra at z = lO/?0, computed 
assuming different values of &o|| Aox = 0.01 (red dashed curve), 0.1 (black 
solid curve), 1 (green dotted curve), 10 (blue dot-dashed curve), and 1000 
(violet dotted curve). 



dissipation. To set the cascade prop ertie s, we utilized the val- 
ues of the constants given in Section lXTl and we also assumed 
s = /i±/a± = 2 and kgn/kox = 0.1. The KAW damping ratio 
j/u> appropriate for the assumed value of (3, which was used 
in Equation dCHK is also shown in green (see also Section|5]l. 
At the outer scale, the peak value of vj_ is 42 km s" 1 . We cau- 
tion that this value should not be assumed to be equivalent to 
the full rms velocity amplitude. In this case, (Ua/ po) 1 ^ 2 = 196 
km s" 1 , which is almost a factor of 5 larger than the maximum 
value of vx at this height. 

The damped spectra shown in Figure [6] have several fea- 
tures that resemble those of measured KAWs in the solar 
wind. Using the conventional form of the reduced energy 
spectrum (eA ~ b\/k±) we found that the magnetic fluctu- 
ation power made a transition from a Kolmogorov-like power 
law k^ 3 to a steeper spectrum with k~^- 5 at k±p p w 1. The 
spectrum becomes shallower again around k±p p rj 40 be- 
cause the wavenumber dependence of (j> flattens out at low 
val ues of j3. This behavi o r is re miniscent of that predicted 
by IVoitenko & De Keysej d201 lh . At larger radial distances 
where j3 > 1, the KAW dispersion relation (Equation d35l )) 
gives rise to a more sustained increase in <f> with increas- 
ing k±. This in turn produces spectra that remain steep, 
with <?a AfJ 2 ' 5 persisting over several orde rs of magnitude of 
k± in agreement wit h both measurements (Smith etal]|2006l; 
Sahraoui et al. 2010) and other models (Howes et al. 2008). 
We note that the predicted undamped KAW power-law de- 
cline of k^ 3 (see Appendix IC.U was not seen for any sus- 
tained range of k± . 

Figure [7] shows the result of varying the normalization 
of the parallel outer scale wavenumber & || on the shape of 
b±(k±). We kept the same value of s = 2 that was used in Fig- 
ure|6] but we varied the constant ratio fcoy / £o_i_ over five orders 
of magnitude. For the lowest values of fe || the outer-scale crit- 
ical balance ratio Xo always remains much smaller than unity. 
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This means that the stirring or forcing takes place well within 
the "filled" region of wavenumber space, and thus strong tur- 
bulence occurs. In this case, b± oc k^ and thus e& oc k^ 3 . 
The opposite extreme case of large & || corresponds to xo ^ 1 
and weak turbulence with less anisotropic driving. In that 

-1/2 

limit, the inertial range spectra are given by b± oc k ± and 
eA k~^. Our model shows the gradual transition between 
these two extreme cases. 

In Figure [8] we compare the Alfven and fast-mode spec- 
tra with one another. As above, we used the background 
conditions at a coronal height of z = 10 Rq and we assumed 
ko\\/ko± = 0.1. We illustrate the most extreme case of a lack 
of high-frequency Alfvenic power by showing the conto urs of 
EpS.k\\,k±) for the case s — > oo. In this limit, Equation dC 16b 
describes an exponential decrease of power with increasing 
X- Other comparable exa mples of this kind of spectrum can 
be found in Figure 4b of ICranmer & van Ballegooijen (2003) 
and Figures 1 and 2 of Uiang et al J (l2009). We computed the 
Alfvenic and fast-mode spectra with the kinetic sources of 
damping that were described in Section [5] Note that £p ex- 
periences the strongest damping at intermediate values of 9. 
For 9 < 10° or 9 > 85°, the transit- time damping described 
by Equation d58l l is relatively weak. 

4. COUPLING BETWEEN ALFVEN AND FAST-MODE WAVES 

4. 1 . Basic Physics and Phenomenological Rates 

There are several ways that the ideal linear MHD wave 
modes can become coupled to one another in the corona and 
solar wind: 

1. Inhomogeneities in the background plasma can blur 
the definitions of the individual modes. For ex- 
ampl e, linear reflection due t o radial variations in 
Va dFerraro & Plumptonl 119581: iHeinemann & Olbertl 
fl980l) " may produce not only incoming Alfven waves 
(i.e., < 1Z < 1), but also fast and slow magnetos onic 
waves (e.g.. lStein|[T97TtlMcDougall & Hoodll2007l) . In 
addition, lar ge-scale bends in the bac kground mag- 
netic field B ilFrischlll964t IWentzelll 19741). density vari- 
ations between flux tubes (IVallevI 11974: Markov!! 
120011: iMecheri & Marschl 120081). or veloc ity shears 
(|Poedts et al.l[l9 98: Gog oberidze et al.l l2007) can drive 
instabilities that partially convert Alfven waves into 
other modes. 

2. Even in a homogeneous medium, the MHD waves 
begin to lose their ideal linear character when their 
amplitudes become large. Nonlinear Alfven waves 
naturally drive second order fluctuations in vii and 
5p that mimic the pr operties of both slow and fast 
magnetos onic wave s (lHollweg|[T97l ISpangleri[l989l 
Vasq uez & Hollwed 119991) . Large-amplitude waves 
also excite a range of wave-wave interactions that 
can often be characterized either as two modes giv- 
ing birth to a third, or one mode sp li tting into several 
others (e.g.. IChin & Wentzej [1971 iGoldsteinl [l978t 
iDel Zanna et alJl200U ISharma & Kumajl2010l) . Mod- 
els of weak turbulence, in which the wave-wave inter- 
actions describe the cascade process (|Chandranl 2005 , 
l2008at iLuo & Melrosell2006t lYooii & FanglfcOOSlT also 
create this kind of coupling. 

3. Although not strictly a multi-mode coupling, when 
k-Lpp ^ 1 the Alfven mode begins to exhibit oscillations 
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FIG. 8. — Comparison of uncoupled power spectra at z = IORq for (a) 
Alfvenic fluctuations, E\(k^ ,kj_), and (b) fast-mode fluctuations, £f(^|| ,k±). 
Contours are plotted one per 10 4 (i.e., one every four decades in power) from 
10~' down to 1(T 29 times the maximum value of £a- Darker shading denotes 
higher power levels. Also shown is a line denoting 9 = 45° (blue dotted curve) 
and the critical balance locus of points that obey Xcff = 1 (red dashed curve). 

in density, parallel electron velocity, and the parallel 
electric and m agnetic fields (Hasegawa & Chenlll976t 
lHollweg| [l999). Observationally, it has proved diffi- 
cult to separate such dispersive KAW density fluctua- 
tions from those arising from independent sources of 
fast or slow MHD wa ves (e.g.. lHarmon & Colesll2005l : 
IChandran et aDl2009h . 

In this paper we take account of one particular nonlin- 
ear ef fect from t he sec ond entry in the above list. Specifi- 
cally, IChandranl (120051) suggested that weak turbulence cou- 
plings between Alfven and fast-mode fluctuations may pro- 
vide enough po wer at high fcn to i nduce substantial ion cy- 
clotron heating. ISuzuki et al.1 (120071) argued that this effect 
may be relatively unimportant because the fast-mode cascade 
timescale Tp is long in comparison to the Alfven cascade 
timescale ta. This may be the case in the low-frequency 
regime of wavenumber space where x ^ 1> but at the cy- 
clotron resonant frequencies of interest (k\\ ~ Q, p /Va) the 
Alfvenic cascade is quenched because x > !■ The fast- 
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mode cascade may in fact even be faster than any intrinsic 
Alfvenic spectral transfer in this region of wavenum ber space. 
Therefore, we proceed using the Chandran (2005 ) results for 
Alfven/fast-mode coupling. 
We express the coupling term in Equations (|3H-(|321 as 

C AF = ^-^ (46) 

7"AF 

such that, in the absence of other processes, the power spec- 
tra at a given wavenumber k are driven toward a common 
value ov er a coupling time scale tafC*). The weak turbulence 
model of Chandran (2005) gave an approximate value for this 
timescale of 

is 

(47) 



7"AF 



15 

23^ 



■tf sin 



which holds in the limiting cases of Ep > Ea and nearly par- 
allel propagation (9 <C 1). In the opposite case of S> Ep, 
it's likely that taf would no longer depend linearly on Tp, and 
may scale instead with r A . However, the region of wavenum- 
ber space with which we are most concerned is the high-fc||, 
low-6* ion cyclotron regime. At those wavenumbers, we know 
that in the absence of coupling the condition Ep ^> Ea is likely 
to be satisfied, and the coupling will be a transfer of energy 
from the dominant fast-mode spectrum to the much less in- 
tense Alfven mode. 
The wave-wave conditions of frequency and wavenum- 



;quency 

ber matching (e.g., Sagd eev & Galeevl [r969) confirm that the 
most rapid coupling should occur when the dispersive prop- 
erties of the Alfven and fast-mode waves are the most similar 
to one another; i.e., at 9 — > 0. Note that Equation (f39b gave 
rp oc 1 / sin 9, so the combined dependence for the coupling 
time is taf oc sin 9. In practice, however, we found that using 
this ideal expression for taf could lead to an unphysical sin- 
gularity at 9 = 0. We removed this singularity by replacing 9 
in Equation ( f4Tb by 9 + 59. We set 69 to a constant value of 
0.01 to avoid having an infinitely fast coupling rate at paral- 
lel propagation. 3 To retain the most generality, we chose to 
reparameterize the coupling timescale as 



r A F = ^t f sin 2 (9 + 6 9) 



(48) 



where we find it useful to vary the constant coupling strength 
$ up or down fr om the value of 237r 2 /15 ~ 15.1 derived by 
IChandranld2005T) . The case <f> = corresponds to ignoring the 
coupling altogether. 

Note that the above form for the coupling timescale implies 
that taf oc k±/k 3 / 2 , so that the coupling is rapid at wavenum- 
bers corresponding to ion cyclotron resonance (large k\\, small 
k±). The coupling is much slower at KAW wavenumbers fa- 
vored by the pure Alfvenic cascade (small kn , large k±). Thus, 
the bulk of the Alfvenic spectrum at x *C 1 is likely to be more 
or less unaffected by the coupling. This seems to be consis- 
tent with our assumption that the integrated energy densities 
Ua and Up also remain uncoupled from one another. We re- 
alize that this may be a severe underestimate of the degree 
of energy transfer between Alfven and magnetosonic modes 
in the corona and solar wind. However, one main purpose of 
this paper is to investigate how much can be accomplished 
with only this small degree of coupling in the high-A:|| tails of 
the power spectra. 

3 Also note that the magnetic field in MHD turbulence undergoes a 
complex, multi-scale "wandering," such that the direction corresponding to 
9 = is continuously va rying in time and space (see, e.g., Ra gotl 120061 ; 
IShalchi & Ko urakis 2007). Thus, the plasma may seldom "see" exactly par- 
allel wavenumber conditions. 



4.2. Approximate Solutions for Coupled Spectra 

The exact solutions to Equations OTb and d32b with cou- 
pling (Caf ¥ 0) must be found numerically. Here we present 
an approximate solution that is both (1) likely to reflect the 
proper behavior of more rigorous numerical solutions in many 
limiting regimes of parameter space, and (2) efficient to im- 
plement on a large grid of model spectra spanning a wide 
range of heliocentric distances. We begin by approaching the 
problem iteratively; i.e., we solve Equation (l3"Tt for Ea under 
the assumption that Ep is known, and we then solve Equa- 
tion ( f32t for Ep under the assumption that Ea is known. The 
analytic solutions derived below suggest a natural way to ter- 
minate this iteration after only one round. 

When solving the advection-diffusion equation for Alfvenic 
fluctuations, let us temporarily ignore the outer-scale source 
term Sa and the dissipation-range damping term that depends 
on 7a- Since we are most concerned with the generation and 
transport of wave power in the high-A;j| regions that undergo 
ion cyclotron resonance, we consider the weak turbulence 
regime of x ^ 1, in which the transport of energy is mainly 
from lo w to high k± and there is negligible parallel spreading 
(see also Oughton et al. 2006). Thus, we solve the advection- 
diffusion equation for discrete, non-interacting "strips" of 
wavenumber space each having constant k\\ . The nonlinear 
coupling supplies wave energy locally, and the Alfvenic cas- 
cade takes it from low to high k±. If we simplify further by 
assuming pure advection (i.e., a± = 0), the time-steady ver- 
sion of Equation OTb becomes 



k i dk i 



TA 



taf 



(49) 



where we use Equation JA6I > to give the timescale ta ~ 
x/(k±y±) in the weak turbulence regime, and we use Equa- 
tion d48l for r A F. 

The above advection-coupling equation can be rewritten as 
a first-order ordinary differential equation, 



dk± 



J0_ + f 

3kj_ 



10/3 



£a = 



f)Ep 

IW3 



where 



/(» = 



$ / v 0F 



\V'0_L 



2 k 1/2 k 5/2 
K 0F K \\ 

,2/ 3 



(50) 



(51) 



To obtain Equation d50b , we made several power-law assump- 
tions for the timescales ta and taf, which depend on the ve- 
locity spectra vj_ (for Alfven waves) and (for fast-mode 
waves), respectively, with 



v_l = V j 



ko± 



-1/3 



Vjfc = V F 



^"0F 



-1/4 



(52) 



We also assumed that we are solving for E A mainly in the 
small-*? region of wavenumber space in which k sa k\\ . 
With the above assumptions taken into account, Equation 
i can be solved by means of an integrating factor. We first 
define the dimensionless independent variable 



y 



3/o 



Ik 



7/3 



7/3 



(53) 



which is a measure of the relative strength of the nonlinear 
coupling. When y ^> 1 (or k± <§; k c ) the coupling is strong 
and we should expect Ea ~ Ep. When y <§; 1 (or k± ^> k c ) the 
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coupling is weak in comparison to the cascade and we expect 
Ea^Ef- Note also that y depends much more sensitively on 
9 than on the magnitude k. Working through the integrating 
factor method and choosing an integration constant of zero 
(to avoid the solution diverging to infinity when y 3> 1), we 
obtain 



45 



ll 

3 



- e y/ ? r 



7' 



(54) 



where T(a,y) is the incomplete gamma function. This func- 
tion behaves as expected in the limits of strong and weak cou- 
pling as discussed above. 

Next we solve the coupled fast-mode advection-diffusion 
equation for Ep under the assumption that E& is known. Mak- 
ing use of many of the same simplifications that were used to 
solve the Ea equation, we include only the cascade and cou- 
pling terms, with 



ap d 



k Dp- 



dEp 
~dk 



Caf 



1 

TAF 



Ep . (55) 



Noticing that the quantity in square brackets above is an ef- 
fective damping rate 7 et f, we use Equation d54l l to write the 
ratio Ea/Ep as a known function of kn and kj_. After substi- 
tuting in the wavenumber dependence for taf, we found that 
7efr oc {k/kop) 1 ! 2 . The analytic solution of Ep(k) for this spe- 
cial case is given in Equation (1C321 >. and the constant c 7 in 
that expression is specified here to be 



' 7 49a F sin 2 6»V" e f 



(56) 



The solution of Equation ( !C32t is applied only for k > kop, 
and the uncoupled/undamped fast-mode power spectrum Eop 
is used for k < kop. 

Since our solution for the ratio Ea/Ep depends only on 
wavenumber and not on any prior solutions of Ea or Ep, we 
found that there is no need for further iteration. We solve 
first for Ep as described above, using Equation (|54l for the 
ratio Ea/Ep, and then we use this ratio to solve for Ea- Note 
that the complete solution for Ep must take account of both 
coupling and transit-time damping (i.e., the damping rate 
given by Equation (l58T>). In practice, we apply both types 
of damping separately to the uncoupled and undamped fast- 
mode power spectrum Eop and we use the solution that gives 
rise to stronger local damping at any given wavenumber. At 
high enough values of feii , the complete solution for Ea must 
take into account the effects of ion cyclotron damping. We 
use the approximate prescription given by Equation ( IC20) to 
implement this damping. 

If the original uncoupled spectra obey Eqa < Eop, then the 
coupled spectra follow 

Eqa < E A < Ep < Eop 

at wavenumbers in the high-fc|| regime where the coupling is 
applied. Usually, the relative increase in Ea from its uncou- 
pled solution is greater than the relative decrease in Ep from 
its uncoupled solution. In all cases, however, we found that 
the variations in the spectra introduced by the coupling do not 
significantly affect the total wavenumber-integrated power in 
either Ea or Ep. 

Figure|9]illustrates the effects of including coupling on Ea- 
As in earlier plots of spectrum results, we used the represen- 
tative height z = IQRq and we assumed £o|| = &o_i_/10. In or- 
der to show that the coupling can be efficient even when the 



40 



w 35 



o 



•t> = l(T 



(a) 



$ = 1 



! * = 10" 



30 



25 



4>=0 



10 



-6 



10 



-4 



10 

k nV n P 



-2 



38 _ 



(b) - 




32 



k„V A /Q p = 10- 



10 6 10 4 10 2 1 10 2 10 4 

FIG. 9. — (a) Slices of time-steady spectra at z = IQRq, shown at con- 
stant k± = kg± : uncoupled spectra Z?oa (black dotted curve) and i?op (red 
dashed curve), and coupled E\ spectra that were computed with a range of 
<E> values (black solid curves), (b) Variation of Ea (black solid curve) and 
£p (red dashed curve) with shown at constant wavenumber k± = an d 
A V\ <->.• 1(1 \ 



uncoupled Alfven wave power Eqa is negligibly small, we as- 
sumed the extreme limiting case of s — > oo. In Figure |9j a) 
we show the kn dependence of the spectra along a slice taken 
at a constant value of k± = ko±. We varied the parameter $ 
between 10~ 6 and 10 +3 . Even if the coup ling is several ord ers 
of magnitude weaker than estimated by Chandra^ (120051) . it 
is still likely to be efficient at generating some Alfvenic wave 
power at kn s» Q p /Va- However, if the coupling constant $ is 
significantly smaller than ^10~ 3 , the ion cyclotron damping 
at m CI p /Va is likely to overwhelm the "local supply" of 
wave energy from the coupling and give rise to a low level of 
resonant wave power. 

Figure |9f b) shows how the power at a given wavenumber 
(k± = ko± and A:y Va/O p = 10~ 3 ) varies as a function of The 
fast-mode power decreases monotonically as $ is increased, 
which confirms our treatment of the coupling in Equation (l55t 
as an effective damping. The Alfvenic power generally in- 
creases (from its uncoupled value far below the lower edge of 
the plot) with increasing <I>, but there is some nonmonotonic- 
ity around $ s» 10~ 2 . This gives rise to a slightly counter- 
intuitive result that there may be more Ea power at high-A:|| 
(and thus more proton and ion heating) at some values of <f> 
than in the $ — > oo limit. 

An example of the full wavenumber dependence of the cou- 
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FIG. 10. — (a) Contours of the £a power spectrum, as in Figure [8]but 
computed with full fast-mode coupling (<E> = 10). (b) Radial dependence of 
spectrum-averaged angle between the background field direction and 
the wavenumber vector k, computed for <E> = 10 (solid curve) and for $ = 
(dotted curve). 

pled £a(£|| , k±) spectrum is shown in Figure [Tol a) for a radial 
distance of r = 10 Rq. This model has the same parameters 
as the one shown in Figure [S] except that we set $ = 10. De- 
spite the appearance of substantial wave power at large values 
of ku, most of the power is still contained within the critical 
balance locus of \ < 1 . This is illustrated in another way by 
Figure [TOl b). in which we show the radial dependence of the 
spectrum-averaged angle Ob* between the background field 
direction and the wavenumber vector k. We used a defini- 
tion for the spec trum-averaged wa vevector anisotropy that is 
similar to that of lGary et alT(l2010h . 



tan 2 (e B *) = 



/ d 3 kE A (k)k\ 
J d 3 kE A (k)kl 



(57) 



Note that the model result at r = 1 AU (89.5°) is reasonably 
close to the value of ^88° measured by Sahraoui et al.l (1201 Ot) 
from the four Cluster satellites at 1 AU. It is evident that 
a strongly perpendicular ("quasi-two-dimensional") sense of 
wavenumber anisotropy is not incompatible with the exis- 
tence of high-frequency ion cyclotron resonant wave power. 

5. KINETIC DISPERSION AND DISSIPATION 

When computing the dissipation rates 7a and tf, we are 
careful to distinguish between two conceptually different 



sources of damping. First, there are the collisional and outer- 
scale cascade processes that were included in Equation (TTTb . 
These processes act at low wavenumber and drive the over- 
all radial evolution of the wave energy densities Ua and Up. 
We do not include them in the damping terms in Equations 
d3Tll-(l32li because their net effects are already included in 
the source terms Sa and Sp. Second, there are the largely 
collisionless kinetic processes that become dominant at large 
wavenumbers. These are the actual processes that dissipate 
the power and give rise to heating, and we describe them in 
the remainder of this section. 

Once the power levels of Alfvenic and fast-mode fluctu- 
ations are specified as detailed functions of ku, k±, and ra- 
dial distance, we compute their damping rates and species- 
dependent heating rates from linear Vlasov theory. Although 
it is known that strong MHD turbulence is far from "wave- 
like" (i.e., coherent wave packets do not survive for more 
than about one period before being shredded by the cas- 
cade), there is a long history of using damped linear wave 
theory to study the small-sc a le dissipation of such a cas- 
cade (s ee, e.g.. | Eichlej 1 1 979t ]Ouataertlll998t iLeamon et al.1 
19991; lOuataert & Gruzinovl 1 1 999t iMarsch & Tul 1200 lal 
Cranmer & van Ballegooijen 2003; iGary & Borovskyl 12004 
Harmon & Colesl 120051) . A typical justification of this ap- 
proach is that no matter the strength of the fluctuations at the 
outer scale, once the cascade reaches the high-A; dissipation 
ran ge the magnitude s are much sma l ler an d quite linear; see 
also lSpanglerfd 19911) and lLehe etafl d2009l) . 

For the Alfven waves, we utilize the Vlas o v-Max well 
code described bv | Cranmer & van Ballegooijen (2002) and 
ICranmer et al.l (120091) to solve the "warm" linear dispersion 
relation for the real and imaginary parts of the frequency in 
the solar wind frame {oj = to r + hf) assuming a known real 
wavevector k. The code uses the Newton-Raphson technique 
to isolate individual solutions from a grid of starting guesses 
in u> r , 7 space, and we select only the left-hand-polarized 
(Alfvenic) solutions. We assumed homogeneous plasma con- 
ditions and isotropic Maxwellian velocity distributions (with 
T p = T e ), and we ran the code for a range of assumed val- 
ues of /3 between 10~ 3 and 10 2 . The code also provides the 
partition fractions of wave energy in electric, magnetic, ki- 
netic, and thermal perturba tions for each wave mode (see also 
iKrauss-Varban et alj|1994h . 

Figure QT| shows several example solutions for the real and 
imaginary parts of the frequency along one-dimensional cuts 
through wavenumber space. For simplicity, we present all 
damping rates 7 as their absolute values, since strictly speak- 
ing the solutions from the Vlasov-Maxwell code all have 
7 < 0. Figure [TTi a) illustrates the approach to the ion cy- 
clotron resonance regime by holding k± constant at a small 
value and plotting u) r and 7 versus ku . Note the cessation of 
weakly damped solutions at 7 w u) r w fl p , which takes place 
at lower values of £|| max for higher values of j3. Equation 
( IC 1 8b is a parameterized fit to the /^-dependence of this cutoff 
wavenumber. 

Figure [TTT b) shows the approach to the high-A:^ KAW dis- 
sipation limit for a constant small value of ku . When solv- 
ing the dispersion relation along a succession of increasing 
values of k±, there are sometimes small discontinuities in 
slope between neighboring solutions (especially in strongly 
damped regions where \j/u> r \ ^ 0.5). Nonetheless, the dis- 
persion properties of our solutions remain sufficiently "KAW- 
like" to represent a continuous set of damping rates from low 
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FIG. 1 1. — Linear dispersion properties of Alfven waves computed for a 
range of plasma /3 values, (a) Real frequencies uj r /Q. p (black curves) and 
damping rates y/u r (red curves) plotted versus hi at constant k±p p = 1CT 3 , 
for = 0.01 (solid curves), /3 = 0. 1 (dashed curves), ft = 1 (dot-dashed 
curves), j3 = 10 (dotted curves), (b) Same quantities as in panel (a), but shown 
as a function of k± at constant V/^/Q p = 10~ 3 . 

to high k±. The behavior of w r versus k± agrees reasonably 
well with the approximate expression given by Equation ( |35V 
For values of /3 > 1, there are secondary maxima in j/uj r at 
k±p P ~ 1 that come from proton Landau damping, whereas 
the larger rates at k±p p > 10 are dominated by electron Lan- 
dau damping. The damping rates shown in Figure [TTT b) were 
also used as the effective KAW ratios ta/^V described in Ap- 
pendix IC.2I These rates were used to compute the high-£ ^ 
dissipation of b± and vj_ as shown in Figures|6]and|7] 

For the fast-mode waves, we make use of a parameterized 
expression for the rate of transit-time damping, which in sev- 
eral studies was found to be the dominant kinetic process to 
dissipate this wave mode (e.g., lBarneslll966t lPerkinslll973l 
lYan & Lazarian 2004). Thus, we assume 



1 + 




exp 



m p /3 cos 2 9 



(58) 



where u r is given by the ideal fast-mode disp ersion relation of 
Equation This expression was given by I Yan & LazariarJ 



(2004) based on initial calculations of Stepanov ( 1958). Equa- 
tion ((58) is valid strictly for only 9 <C 1 , but it does not diverge 
from the more exact solution at larger 9 by more than about a 
factor of two. 

The remainder of this section describes how the dissipated 
Alfven wave energy is partitioned between protons, electrons, 
and heavy ions. We ignore the particle heating that comes 
from fast-mode wave dissipation because its overall magni- 
tude was found to be small in comparison to that from Alfven 
waves. In a pure hydrogen plasma, we separate the damping 
rate 7 into components attributed to the kinetic effects of pro- 
tons and electrons. To zeroth order, the contribution to 7 from 
other ions is negligibly small and can be estimated separately 
(see below). Thus, we define 7 = 7^ + 7^, where 



7i = 7 



(59) 



where s = p,e denotes either the protons or electrons, and the 
species-dependent resonance functions are given by 



2 +00 
i>s = 1— > , eX P 



Fll« 



£=-00 



m=e-\ 



l mhn{(,±) 7 



(60) 



where w 2 s = Ane 2 n s /m s is the squared plasma frequency, 
vviu and w± s are parallel and perpendicular thermal speeds 
of species s, and /,„ is the m-order modified Bessel func- 
tion of the first kind. The dimensionless coefficients 
a m depend on the electric-field polarization vector that 
is output from the Vlasov - Maxw ell dispersion code of 
ICranmer & van Ballegooiien (200 3]), and th e y are g iven in 
full by Equations (43)-(45) of iMarsch & Tul feOOlah . Equa- 
tion ( f60b is valid for an isotropic Maxwellian distribution, for 
which W\\ s = w± s and there is assumed to be zero differential 
bulk flow between the protons and electrons. The dominance 
of ion cyclotron or Landau damping depends on the values of 
the dimensionless resonance factors, 



■m. 



k\\w\\ 



k±w ±s 



(61) 



In practice, we truncate the infinite sum in Equation (1601 at 
— 10 < £ < +10. Test runs made with a larger range of sum- 
mation indices produced no substantial differences from those 
using the default range. 

Figure [L2l shows separate sets of contours for 7 P /w r and 
je/u r in wavenumber space for an example value of (3 = 
0.1. These contours can be comp ared with Figure 4(a) of 
ICranmer & van Ballegooiien (2003). which was computed for 
f3 pa 0.01. The proton damping rate jp/u,- increases rapidly 
as k\\VA/Q p approaches unity, and the electron damping rate 
je/uj,- increases more slowly as k±p p increases from 0.1 to 
100. The complex behavior of the contours in region of 
wavenumber space with both high k\\ and high k± is the re- 
sult of the dispersion relation being affected by th e presence 
of strongly dampe d ion Bernstein modes (see, e.g.. lStrxlll992l : 
iHowes et al.ll2008l) . 

It is evident from FigureQ~2]that, in the solar corona, the re- 
gion of nearly parallel Alfven wave propagation in wavenum- 
ber space (i.e., 9 <C 1) is dominated by proton damping and 
the region of nearly perpendicular propagation (9 — > ir/2) is 
dominated by electron damping. The o bservational ev idence 
for preferential proton and ion heating dKohl et al.ll2006l) thus 
presents a problem when confronted with the dominant per- 
pendicular anisotropy of Alfvenic turbulence. 
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FIG. 12. — Contours of ■y p /co r (thick curves) and y e /to r (thin curves sepa- 
rated by varying gray shading) plotted versus and kj_ . Contours are plotted 

twice per decade from 3 X 10 -5 to 3 X 10~' and generally go from low to high 
values with increasing wavenumber. A line denoting 8 = 45° (dotted curve) 
and a point illustrating where 8 ep is defined (filled circle) are also shown. 
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FIG. 13. — Radial dependence of the tangents of 8 ep (solid purple curve) 
and 6> cr j t (black curves), the latter computed for Ar || /^o_L = 01 (dotted), 
^o]|Ao_L =0.1 (dot-dashed), and £ || Ao_L = 1 (dashed). The gray region de- 
notes the approximate region of parameter space expected to be "occupied" 
by a purely Alfvenic turbulent cascade. 



Figure[T3]illustrates the magnitude of this apparent discrep- 
ancy by comparing the large-scale radial dependence of two 
key angles. The strongly anisotropic Alfvenic cascade is illus- 
trated by 9 a i t , w hich is the angle between k and Bo at which 
occurs both the Goldreich & Sridhar (1995) critical balance 
(X = 1) and the onset of KAW dispersion (k±p p = 1). We find 



that tanf? cl -it ~ Va/^±, where b± is evaluated at k±p p = 1, and 
we plot tan^crit for three example values of the outer-scale 
wavenumber ratio k ^/ko±. Figure [13] also shows the radial 
dependence of 9 ep , which is defined as the angle at which the 
contours for 7 P /w = 0.1 intersect with those of 7 e /w = 0.1 in 
wavenumber space. (This point is shown in Figure [T2] with a 
filled circle.) For 9 < 6 ep the damping is dominated by protons 
and ions; for 9 > 9 ep the damping is dominated by electrons. 
Note that 9 ep <C # C nt in the solar corona and much of the inner 
heliosphere, so that it is difficult to see how the cascade of 
linear Alfven waves alone can be responsible for the observed 
proton and ion heating. 

We computed the rates of proton and electron plasma heat- 
ing from the modeled values of 7 P and j e by usin g the 
quasilinear framework outlined by | Marsch & Tul (1200 lab and 
ICranmer & van Ballegooijen (2003). The volumetric heating 
rates Q s (e.g., expressed in units of erg s" 1 cm" 3 ) are given by 
integrals over vector wavenumber k of the form 



9 



— = / £/ J k£ A (k)2 7 , 



(62) 



where s = p,e denotes the particle type of interest. For 
now, we ignore differences between parallel and perpendic- 
ular heating and only compute the summed heating rate Q s = 
Qs\\ +Qs±- In order to perform the wavenumber integration 
in Equation d62l l. we constructed two-dimensional numerical 
grids of 7 P and j e for values of 10~ 3 Q p /Va < &|| < £|| max and 
10" 3 < kxPp < 10 3 . We used 200 points in k\\ and 100 points 
in k±, and we constructed a total of 14 grids for values of /3 
ranging from 10~ 3 to 22 (with f3 varying logarithmically with 
three samples per decade). Linear interpolation was used to 
evaluate the damping rates at values of kn, k±, and (3 between 
the discrete grid points. We assumed that the ratios j p /oJ r and 
%/u) r remain constant as one extrapolates into the weakly- 
damped regions defined by k±p p < 10~ 3 and A;||VA/f2p < 10~ 3 . 

To estimate the heating rates experienced by heavy ions, 
we assume that most low-abundance ions do not have a sig- 
nificant effect on the overall wave dispersion relation. This 
allows us to use an "optically thin" resona nce condition fo r 
the ion cyclotron wave-particle interaction (ICranmerl 12000), 
which results in a perpendicular heating rate 

Qn _ Trfif Zi\ f d 3 kEA(k)S(k ^/Va), (63) 



V A 



ion charg e and mass i n pro- 
ra nnierl 1200 11 iTu & Marsclil [2001; 



where Z, and A, ar e the 

ton units (see also | Cran ._ ... .... ... _ 

lLandi & Cranmerl 12009). The Dirac delta function extracts a 
one-dimensional "strip" of the power spectrum that is in res- 
onance with the ion Larmor motions at uj r « feii Va = Thus, 
Equation d63l can be evaluated with just a single integration 
along the k± direction. 

6. RESULTS FOR COLLISIONLESS PARTICLE HEATING 

Here we present results for Q p /Q e , the ratio of proton to 
electron heating rates, computed from Equation (l62l with var- 
ious assumptions for the shape of the turbulent Alfven-wave 
spectrum E A (k«,k±). Figures [141 and [TBI show how this ra- 
tio behaves for pure "uncoupled" Alfven waves, and Figure 
[16] summarizes the outcome of coupling the Alfven and fast- 
mode waves as discussed in Section [4] Table [2] summarizes 
the specific values of cascade and coupling parameters that 
were assumed in each of these plots. 

In Figure [14] we show the radial dependence of Q p /Q e 
for various methods of computing the outer-scale parallel 
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FIG. 14. — Radial dependence of logQ p /Q e for: (a) constant Alfven wave 
periods P = 1 min (solid red curve), P = 10 min (green dotted curve), and P = 
100 min (blue dot-dashed curve); (b) outer-scale & || determined from ideal 
Goldreich & Sridhar ( 1995) critical balance (red solid curve) and from mod- 
ified Beresnvak & Lazarian (2008, 2009) critical balance (blue dot-dashed 
curve); (c) constant ratios k ^/koj_ = 0.01 (red solid curve), ^o|| Ao± = 0.1 
(green dotted curve), &6||/&ox = 1 (blue dot-dashed cu rve), and k ^/kp ±^ = 10 
(purple dashed curve). Also shown in (a)-( c ) are the ICranmer et aTJ ( 2009) 
measurements (gray region) and the Howes (2010) model prediction (black 
dashed curve). 



wavenumber fe ||' In all panels, the turbulence spectra were 
computed with constant values of s = 2 and $ = 0, as well 
as the other default parameter choices discussed in Section 
13.11 Figure PHI" a) assumes a range of radially constant wave 
frequencies which d etermine fco || from Equation (1431 . Fig- 
ure [l4|b) applies the Goldreich & Sridhar (1995) conditions 
of critical balance for both zero and nonzero cross helicity at 
the outer scale; see Equations (|441>-(|45T>. 

Figure PBl c) shows the relative heating rates Q p /Q e for 
a range of constant ratios k ^/ko±. At the coronal base 
(z = 0.01 Rq), note that Q p /Q e behaves non-monotonically as 
a function of this wavenumber anisotropy ratio. The mini- 
mum value of Qp/Qe occurs at &o||/&oj_ ~ 0.55. The non- 
monotonic behavior occurs because of two competing effects. 
At large values of fc ||> the weak-turbulence critical balance 
curve (xo = 1) begins to approach the ion cyclotron frequen- 
cies. This has the result of increasing Q p while leaving Q e un- 
changed. However, when fe || becomes very small, the wave 
power becomes concentrated into narrower "cones" that pro- 
vide more energy to the KAWs. This has the result of increas- 
ing both Q p and Q e , but the smaller rate Q p receives a larger 
fractional change. 

Each panel of Figure [14] also shows the empirically deter- 
mined range of Q p /Q e rati os from the Helios a nd Ulysses 
measurements described by Cran meret al.l (12009b . The plot- 
ted error range of ±0.3 in \og{Q p /Q e ) accounts for both 
modeling and observational uncertainties. Also, we show 
the theoret ical predictio n for Q p /Q e from the gyrokinetic 
model of Howes (2010) as a dashed black curve. As dis- 
cussed by Howes (2011), this model agrees well with the 
ICranmer et al.l ([2009) measurements at r > 200 .Rq, but un- 
deresti mates the proton heating at r < IQQRq. The iHowesI 
(2010) gyrokinetic model includes the same sources of high- 
Ic± KAW damping that we use, but not the high-A:|| sources 
of ion cyclotron dampin g. In Figure [I4l we find that the 
best agreement with the lCranmer et al.l (12009) measured ratio 
comes from the model that assume s critical balance with the 
iBeresnyak & Lazarianl (12008. 2009) modification for nonzero 
cross helicity; i.e., xo ~ l/7£. 

In Figure [15] we vary the ratio s used i n the Alfvenic paral- 
lel cascade function g(x)', see Equation ( IC14b . We retain the 
Xo ~ 1 /H approximation for k^n that was found to be an opti- 
mal choice for agreement with observations at r > 6QRq. For 
lower heights in the low-/? corona, we find that large values 
of s give insufficient wave power at the ion cyclotron reso- 
nant values of k\\ to provide significant energy to the protons. 
One would need to specify s < 0.5 in order for there to be 
enough high-A: | power t o give protons a substantial f r action 
of t he dissipated energy. ICranmer & van Ballego"o iien (2003) 
and lLandi & Cranmerl (120091) came to this same essential con- 
clusion. Although there are still no firm experimental or theo- 
retical bounds on the expected value of s in MHD turbulence, 
it is generally believed that values as low as s < 0.5 are unre- 
alistic. 

Figure [16] shows the results of mode coupling between the 
Alfven and fast-mode fluctuations. The curves in Figure \l6\ a) 
were computed for a range of constant values of the coupling 
constant $ from 10~ 6 to 10 +3 . At large distances (r > 0.3 
AU), it is clear that the presence or absence of coupling has 
very little effect on the Q p /Q e ratio. This insensitivity occurs 
because much of the proton heating at intermediate and high 
values of j3 comes from the Landau and transit-time damping 
of KAWs. The low-kn, high-A:^ part of the spectrum is 
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FIG. 1 5 . — Radial dependence of log Q p /Q e for xo ~ 1 /H and s = 0-25 (red 
solid curve), s = 0.5 (orange dashed curve), s = 1 (green dot-dashed curve), 
s = 2 (cyan solid curve), s = 4 ( dark blue dotted curv e), and s = 8 (black 
dashed curve). Also shown are the Cranmer et al. (2009) measurements (gray 
region). 

there no matter the value of and it dominates the proton and 
el ectron heating in this case. The results are similar to those 
of lHowes] d2010ll20Tlh who did not include mode coupling. 

In the low-/3 corona, Figure [T6l indicates that $ needs to 
be at least of order unity to excite sufficient power in high-£|| 
ion cyclotron waves to heat protons on par with the electrons 
(i.e., Qp/Qe ~ 1)- For low values of $, the plotted ratio un- 
dergoes several increases and decreases as a function of radius 
that we cannot trace to any one simple cause. The local maxi- 
mum that appears at z ~ 0.5 Rq corresponds to the local min- 
imum in plasma j3 (see Figure [TJ. In the low-/? regime, it is 
likely that the relative "competition" between mode coupling, 
transit-time damping (for Ep), and ion cyclotron damping (for 
£a) undergoes numerous reversals as a function of radius. 

In Figure [ToT b) we fix the coupling constant at $ = 10, 
which is of the same order of magnitude as predicted by 
IChandranl d2005l) . and we vary the normalization of the fast- 
mode wave power. It was evident from Figure |4] that small 
changes in the large-scale wave transport properties could 
give rise to large changes in the fast-mode power in much 
of the corona and solar wind. Thus, we take the standard 
model for Up(r) and multiply it by constant factors ranging 
from 10~ 3 to 10 +2 . We note, however, that we do not have 
excessive freedom to increase the Up normalization too far 
above the standard model. A significantly higher coronal pop- 
ulation of fast-mode waves would contribute to a larger vj_ 
that may exceed the observational constraints shown in Fig- 
ure |3l a). Nonetheless, Figure [ToT b) shows that the standard 
model ends up being a reason able solution that m atches the 
observed in situ heating ratio (Cra nmer et al.ll2009b and also 
gives appreciable proton heating in the extended corona (as 
required qualita tively from UVCS proton tempera ture mea- 
surements); see ICranmer & van Ballegooiien (2003|). 

An example calculation of preferential heavy ion heating 
is shown in Figure [T7] The ion used for the model was 
+5 , whose properties have been measured in the corona 
from emission in the O VI 1032, 1037 A spectral line doublet 
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FIG. 16. — Radial dependence of logQ p /Q e for varying properties of 
Alfven/fast mode coupling, with: (a) standard model for 1/f and a range 
of coupling constants: $ = (red solid curve), <3? = 10~* (orange dashed 
curve), <!> = 10~ 3 (green dot-dashed curve), <E> = 1 (black solid curve), <E> = 10 3 
(dark blue dotted curve); (b) constant value of <J> = 10 and a range of mod- 
ified values for fast-mode power: t/ F /10 3 (dark blue dotted curve), U r /10 2 
(cyan solid curve), t/p/10 (green dot-dashed curve), the standard model of 
J7f (black solid curve), 101/f (ora nge dashed curve), 100 [/f (solid red curve). 
Also shown in both panels are the Cranmer et al. (2009) measurements (gray 
regions). 

dKohl et al]|2006l) . We used the parameters corresponding to 
the best agreement with observational constraints on Q p /Q e 
(see Table |2|. We then adjusted the fast-mode wave power 
Up(r) by changing the multiplicative constant that was varied 
in Figure [ToT b). As in Figure [ToT b), values of this multiplica- 
tive constant between about 1 and 10 appear to bracket the 
observational constraints. 

The plotted ranges for the observationally determined Qxi 
rates w ere derived by combining observatio ns from both the 
UVC S ( ICranmer et al.fl999l) and SUMER dLandi & CranmerJ 
2009) instruments on SOHO with semi-empirical solutions 
of the perpendicular internal energy conservation equations. 
These heating rates we re not given exp l icitly by either 
ICranmer et ail (1 1999b or lLandi & Cranmer! (120091) . but they 
were computed and saved from the models that produced 
agreement with the observed radial behavior of 7jj. The 
SUMER and UVCS data were obtained for off-limb measure- 
ments of O VI emission, in which the line widths are primary 
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FIG. 17. — Radial dependence of the perpendicular heating rate per unit 
mass Q±i/(ntirii), in units of erg s -1 g" 1 , for +s ions. Model results shown 
for a range of modified values for fast-mode power: J/f/100 (dark blue dot- 
ted curve), [/f/10 (green dot-dashed curve), standard Uf (black solid curve), 
IOCf (orange dashed curve), 100C/f (solid red curve). Also shown are empir- 
ical constraints from SUMER and UVCS emission line measurements (gray 
regions). 

diagnostics of 7jj. Note that the radial dependence of the 
two observationally determined regions is similar to that in 
the plotted model curves. However, the SUMER data corre- 
sponds to about a factor of 10 higher fast-mode wave power 
normalization than the UVCS data. 

If the postulated mode-coupling explanation for ion cy- 
clotron proton/ion heating is correct, then the results given 
in Figures [16] and [17] constrain the required levels of fast- 
mode wave power. In the low corona (z < 0.1 R®), there may 
need to be up to a factor of 10 higher value of U-p than in the 
standard model of Section [2] but in the extended corona and 
heliosphere the standard model may be close to correct. Of 
course, it is only the high-£|| tail of the fast-mode spectrum 
that matters to the calculation of available Alfvenic power at 
the ion cyclotron resonances, not its outer-scale normaliza- 
tion. Therefore it is possible that Up may depart significantly 
from the values predicted by the standard model of Section 
[2] but still produce agreement with the various observations 
by having different values for the spectral slope and angle- 
dependence of E-p(\i.). 

7. DISCUSSION AND CONCLUSIONS 

The aim o f this paper was to explore the consequences of 
[Chandrae's (12005b conjecture that nonlinear couplings be- 
tween Alfven and fast-mode waves may produce sufficient 
ion cyclotron wave power to heat protons and heavy ions in 
the corona. To test this idea, we constructed a semi-empirical 
model of the background plasma and MHD wave properties in 
a flux tube connected to a polar coronal hole. For the sake of 
practicality, we utilized several approximations when solving 
the wave energy transport equations for the energy densities 
of Alfven, fast, and slow modes: 

1 . The equations themselves were adapted from standard 
WKB "wave action conservation" theory, which does 
not take account of the effects of linear wave reflection 



in a fully self-consistent manner. We also assumed the 
associated WKB limiting case of equipartition between 
the kinetic and magnetic energy densities for the dom- 
inant Alfven waves (i.e., K y = M r ). Roughly speaking, 
these approximations are consistent with an assumption 
that the wave frequencies are higher than ~ 10~ 3 Hz 
in the corona. However, it has also been shown that 
the radial behavior of Alfvenic wave power in the so- 
lar wind is never far from the predictions of WKB the- 
ory eve n in the heliosphere where reflection is not neg- 
ligibl e (IZank et al.lll996t ICranmer & van Ba llegooiien 
l2005h . 

2. Because of other evidence that the dominant inertial- 
frame frequencies in coronal MH D turbulence may be 
lower than ~ 10~ 4 H z (see, e.g.. IChandran & Hollwegl 
2009: Cranmerll20l0l) . we made use of a low-frequency 
approximation for the Alfven wave reflection coeffi- 
cient 1Z. This also involved an analytic approxima- 
tion for the radial dependence of the Alfven speed scale 
height Ha (Equation (f25ll). 

3. For the fast and slow magnetosonic waves, we modeled 
the radial transport of an isotropic ensemble of propa- 
gation directions 6 using a single wave action conserva- 
tion equation. We chose one reasonable method to per- 
form the averages over 9, but other methods may yield 
different results. We also used the Eulerian average for 
the outflow speed uq and neglected the second-order 
effects of "Stokes drift" that would enter int o the as- 
sociate d Lagrangian version of the mean (see Cranmer 
l2009bh . 

Although the effects of removing these approximations 
should be investigated further, we do not believe their use in- 
validates the results of the wave transport models presented 
above. 

With the above caveats taken into account, we produced a 
standard model of the Alfven, fast, and slow mode energy 
densities between 0.01 and 1000 Rq above the solar photo- 
sphere. In agreement with earlier results, we found that slow- 
mode MHD waves of solar origin probably cannot survive 
into the extended corona. In addition, we found that the am- 
plitudes of fast-mode waves at large distances are more sen- 
sitive to the assumed model parameters than are the ampli- 
tudes of Alfven waves. For this reason the standard model of 
fast-mode wave energy density was treated as a representative 
example and not a definitive prediction. Thus, other reason- 
able models of the available fast-mode power can be obtained 
by multiplying or dividing the standard model's energy den- 
sity by factors of order 10-100 without sacrificing too much 
realism. 

At each radial distance, we simulated the time-steady 
wavenumber power spectra of Alfvenic and fast-mode turbu- 
lent fluctuations. We included the effects of nonlinear cou- 
pling and collisionless kinetic wave dissipation. We also 
computed the time-steady heating rates for protons, electrons, 
and a representative minor ion species (0 +5 ) for comparison 
with observational constraints. The resulting heating rates for 
the standard model of fast-mode wave power was found to 
provide both substantial heating for coronal protons as well 
as produce agreement with the preferential O" 1 " 5 ion heating 
measured by WCS/SOHO. However, if the fast-mode wave 
power in the corona is significantly lower than was assumed 
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in the standard model, the proposed idea of mode coupling is 
probably not a viable mechanism for the ion heating. 

In order to match some of the observations — such as the 
need for Q p /Q e to be of order unity at z < 0.1/?© and 
for the +5 heating rate to agree with that measured by 
SUMER/SOHO at similar heights — we found that approx- 
imately 10 times the standard model's assumed fast-mode 
wave energy density may need to be present. This could 
be accounted for in several ways. First, we neglected the 
effects of Alfven waves giving rise to second order fluctua- 
tions that mimic the properties of both fast and slow mag - 
netosonic waves dHoliwedll97lt IVasauez & H6iiwegl [T999). 
It is possible that these secondary oscillations could behave 
similarly enough to ideal fast-mode waves that they enable 
the same kinds of cascade and coupling. Second, we also 
neglected nonlinear couplings that involve slow-mode MHD 
waves, which appear to dominate the density fluctuations in 
the l ow corona. It may be possible for these couplings (see, 
e.g.. lYoon~& Fang 2008) to also power the high-A:|| part of the 
Alfvenic fluctuation spectrum. 

To make further progress with the proposed set of ideas, 
it will be important to better understand the origin of the 
fast, slow, and Alfven w a ves in t he solar photo sphere and 
chromosphere. iHollwegl (Il978al) . ISpruitl (1 19811) . and oth- 
ers studied the wavelike oscillations induced by convective 
jostling in small-scale flux tubes that extend up into the 
chromosphere. However, once waves reach the sharp and 
"corrugated" TR boundary, they can undergo refle ction, re- 
fraction, and multiple types of mode conversion dHp llweg 
1978bt iBogdan et all 120021; lHasan & van Bal legooiienl l2008t 
Fedun et al] 1201 U ICallv & Hansen! 1201 lb . The types and 
strengths of MHD waves that survive the chaotic lower at- 
mosphere probably also depend on the nature of the region 
underlying the solar wind flux tubes of interest (i.e., coronal 
hole, active region, or quiet loops). 

Future work must also involve more physical realism for 
the model of turbulent cascade. Replacing our hodge-podge 
collection of analytic solutions with a fully self-consistent 
numerical simulation is an obvious priority. A key part 
of this improvement will be to remove the assumption of 
scale separation that prevents different radial zones from in- 
teracting with one another in wavenumber space (see, e.g., 
IVerdini et al.ll2009t) . In addition, we note that the advection- 
diffusion terms in Equations d3~TT>-(l3"2li contain the limiting 
assumption that the spectral transfer is "local" in k-space. It 
has been shown that true MHD turbulence is not so local be- 
cause of intermittent high-order wave- wave interactions and 
nonlinear steepening effects (e .g., Medvedev 2000; Chandran 
2008b; ICholl2010t iHowes et alj|2012h . We also assumed en- 
ergy equipartition between the vj_ and b±_ spectra in the MHD 



inertial range, but in situ measurements show that not to be 
the case in actual solar wind turbulence (Grappin et al1 ll983t 
IWang et alJ20H . 

We also intend to improve upon the kinetic treatment of col- 
lisionless particle heating described in Section|5] We assumed 
isotropic Maxwellian velocit y distributions whe n solving for 
the linear damping rates, but Bash ir et ail (1201 Ol) showed how 
non-Maxwellian temperature anisotropies can significantly 
affect the KAW dispersion relation. The ultimate rate of elec- 
tron heating from KAW Landau damping can also be affected 
by nonlinearity a nd Coulomb collision e ffects that we did 
not include (e.g.. iBorovsky & Garyll20TTl) . The time evolu- 
tion of proton and ion velocity distributions, under the in- 
fluence of cy clotron resonant heating, is also decidedly non- 
Maxwellian dGalinskv & Shevchenkoll2q00t llsenberd l200Tt 
ICranmerll200lHlsenberg & Vasquezll2009ll201 ll) . 

Finally, we emphasize that the proposed idea of nonlin- 
ear coupling between Alfven and fast-mode waves is only 
one proposed solution to the problem of preferential pro- 
ton/ion heating. Some of the other suggested explanations 
were listed briefly in Section [TJ One recent example that has 
received significant attention is the stochastic energization of 
protons and ions t hat occurs when KAW amplitudes become 
sufficiently high dJohnson & Chengj l200lt IChandranl l2010t 
IChandran et al j|201 ll) . To excite this proposed stochasticity, 
the dimensionless ratio v±/c s (evaluated at kj_p p =1) should 
exceed values of order 0. 1 . However, in this paper's standard 
model of Alfvenic fluctuations (either with or without non- 
linear couplings), this ratio never exceeds a value of 0.003. 
The main factor re sponsible for this dramatic mismatch is our 
assumption of the iGoldreich & Srid har (1995) scaling in the 

-1/3 

limit of strong turbulence (i.e., v_l cx k ). Alternate theories 
of the strong Alfvenic cascade (e.g., Boldvrev 2006t IPodestal 

give a shallower dependence of vj_ cx kj_ . This would 
allow larger values of vj_ to survive to the onset of KAW dis- 
persion at k±p p w 1. We await improved theoretical descrip- 
tions of MHD turbulence and conclusive empirical tests of 
such models. 
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APPENDIX 

A. HEURISTIC OVERVIEW OF MHD TURBULENCE 

The cascade of energy from large to small eddies was first described in the co ntext of isotropic hydrodynamic turbulence 
dvon Karman & How arth 1938; Kolmo gorovll 1 94 ll ; lObukhovl 1 94 lHBatchelor| [l95 3 ) . The spectral transport timescale for energy 
to be transferred down to the next order of magnitude of eddy size is estimated generally as t s {kvkT 1 , where k is the magnitude 
of the local wavevector k and is the local eddy velocity at this value of k. For isotropic fluctuations that depend only on k and 
not its direction, we can define the reduced one-dimensional spectrum e m (k) = v\/k. Thus, since 



Po 



— / die 



(Al) 



we relate the eddy velocity to the full three-dimensional spectrum via v\ = \-nk?E m . The cascade rate is estimated as e ~ v\j 
Assuming that e is constant in the inertial range leads to the time-steady Kolmogorov-Obukhov spectrum e m cx £~ 5 / 3 , or E m 
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r 1 '/ 3 . 

When the bac kground magnetic field becomes strong, other physical processes become important. Iroshnikov (1 19631) and 
iKraichnanl d 19651) (hereafter IK) realized that the "eddy" description of hydrodynamic turbulence could be generalized by re- 
ferring to colliding MHD wave packets, and that the Alfven speed Va introduces a new absolute scale into the problem. If one 
continues to treat the cascade isotropically in £-space, a more generalized spectral transport time can be defined as 

(A2) 



kv k \ v k 

where p = gives the Kolmogorov-Obukhov limit and p = 1 is the result of the IK analysis. Using the same assumption above 
that e is constant, we obtain a more general one-dimensional power spectrum e m oc k~ { - p+5) l {p ^ ) . For the IK value of p = 1, the 
spectrum is <?,„ oc A;" 3 / 2 (see also lBoldvrevl2005h . 

It has been known for several decades that a cascade of Alfven-wave-like fluctuations does not lead to an isotro pic distribution 
of power in wavenumber space (IStraussll 19761 iMontgomerv & Turner|[l98lt IShebalin"eFaI|[T98l lHigdonlfT984h . The dominant 
energy cascade takes place mainly in the two-dimensional plane perpendicular to the background field. For the Alfvenic fluctu- 
ations, we can define the local eddy velocity as v_i_ being mainly a function of k±. The one-dimensional spectrum in this case is 
given by e& = v 2 L /kj_ and the integration over wavenumber space is best done in cylindrical coordinates with 

U A f f f 

= / dkn I dk±k±EA = / dk±e& ■ (A3) 

271-po J J J 

Taking into account the spectral anisotropy (kn =/k±)we can also write an even more general perpendicular transport time for the 
Alfven waves as 

A k±v± (vj_) (&J_) *" 

A perpendicular generalization of th e IK model is given by p = 1 and q = 0, which gives e A oc k^ 2 (see also lNakavamalfT999l 
2001; Boldyrev 2006; Podesta 20 l ib. Weak three -wave coupl i ngs have been shown to give rise to the case p = q= 1, which yields 
e\ oc k~^ (e.g.. iGaltier et al.l l2000; Bhattach ariee & Ngll2001t iBoldvrev & Pere z 2009). However, in that case nonlinear effects 
grow in magnitude as k± gets larger, so it is generally believed that a weakly turbulent inertial range must eventually become 
st rongly turbulent (see a lso |Goldreich & Sridharll 1997b . 
IGoldreich & S ridhar ( 19951 described strong Alfvenic turbulence with a spectral transfer time given by p = q = 0, and thus 

-5/3 

ca oc k ± reminiscent of the original Kolmogorov-Obukhov model. In this case of strong mixing between the turbulent motions 
(perpendicular to the field) and the flow of Alfven wave packets (parallel to the field) there is a so-called "critical balance" that 
couples k± and k\\ to one another. We define a critical balance parameter 

X « f-=- (A5) 

k±_v±_ 

such that the IGoldreich & Sridharl (fl995) strong cascade is consistent with the condition \ s» 1. Combining this with the velocity 
scaling v± oc k^ yields the wavenumber anisotropy scaling fc|| oc k 2 / 3 . Note that assuming p = q in Equation ( IA4b is equiva- 
lent to ta being given by \ p /(fej_vj_); see also IGaltier et al.l d2005). An alternate way of describing the cascade was given by 



_JLX . 

IZhou & Matthaeusl dl990bl) . who defined a triple correlation timescale equivalent to 

1+v 

ta « ~. . (A6) 

k±_v±_ 

This ex pres sion naturally bridges the strong (x < 1) and weak (x ^> 1) turbulence scaling limits, and we use a similar relation in 
SectionO 

The cascade of compressible fast-mode waves has received less attention than that of the incompressible Alfven waves. Because 
fast-mode waves propagate at a roughly constant phase speed no matter the direction angle 8, the fast-mode cascade has been 
suspected to resemble isotropic hydrodynamic turbulenc e. In fact, numerical simulations do tend to fi nd that fast-mode waves 
produce a more isotropic spectrum than do Alfven waves (ICho & Laz arian 2003t lSvidzinski et alj2009|). 4 The rate of the cascade 
is gen erally assumed to follow the weak IK-type scaling of Equation (1A21 > with p= 1 (see. e.g.. lChandranll2005l; ISuzuki et al.l 
2007). Thus, because in most cases we expect <C Va. the fast-mode cascade timescale rp is likely to be significantly longer 
than the Alfvenic timescale ta. 

It is important to emphasize that there is still no agreement concerning the most realistic and universal way to describe MHD 
turbulence. There remains cont roversy about the applicability of the various power law exponents (es pecially 5/3 versus 3/2) 
for the Alfvenic inertial range (Beresnyak 2011; Forman et alJl201 it iMason et alJl20TH : lPodestall201 ll) . Simulations have not 
been able to accurately pin down the amount of slow "leakage" of power to the high-A:|| region of the spe ctrum where y 1 . 
Furthermore, the observed steepening o f the spectrum at high values of k± is still not well understood dLeamon et al.lll998t 
IStawicki et afll200H: IMowes et alJl2008h . In many models, the precise scalings depend on the degree of cross helicity of the 
fluctuations (i.e.. on the imbalance betwe en Z+ and Z_) and on whether the turbulence is driven or decaying (e.g.. iLithwick et al] 
120071; IChan dran 2008b; Chen et aTll201 ll) . In this paper, we attempt to identify the most controversial aspects of the models and 
discuss how they can be modified once such issues are resolved. 

4 This is generally valid in the ideal MHD range, at which ui < f2 p . We ignore the large literature of "whistler turbulence," in which dispersive effects may 
lead to wavenumber anisotropy at higher frequencies ui ^> 
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B. LINEAR DAMPING RATES: COLLISIONAL AND COLLISIONLESS 

lAlfvenl d 19471) and[Osterbrock| (119611) first proposed that MHD waves in the solar atmosphere could be damped by collisional 
processes. These processes include viscosity, thermal conductivity, electrical resistivity (i.e., Joule or Ohmic dissipa tion), and 
ion-neutral friction. In the fully collisional regime we make use of the basic expressions derived by Braginskii (1965). Here we 
describe the total linear damping rate for MHD wave mode m by a sum of three components, 

7'« — 7vis,m + 7ohm,m + 7con,m (Bl) 

where 7 v i s .m denotes damping due to kinematic viscosity, 7 hm,m denotes electrical resistivity, and 7 con .„, denotes thermal conduc- 
tivity. Since our main goal is to model the wave damping in the (almost completely ionized) corona and solar wind, we ignore 
ion-neutral friction. For Alfven waves, 



7vis,A 
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— (m*±+»fc*jO (B2) 
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4tt<t\\ J \ Attct j 



7ohm,A = { t— I *i+ - A k\ (B3) 



7con,A = . (B4) 

For fast and slow mode waves (m = F,S), we note that the damping rates given explicitly by Braginskii ( 1965) are valid only in the 
/3 <C 1 limit. The expressions given here are appropriate for arbitrary values of f3, but we made the assumption that k\\ ~ k±. In 
other words, for the assumed isotropic distribution of fast and slow wave vectors, the damping rates depend only on the magnitude 
k 2 = kl + k 2 L . With that caveat, the damping rates are given by 



2k 2 r % 



Po 



y (4/« - 4^ + /,, v ) + mf vx + m (/„ + 2f xz + f vz ) 
c 2 k 2 f B 



(B5) 



7ohm,m = -j (B6) 

47TO"j_ 

4( 7 -i)^ B r 2( . 

5 fthk {K\\+K±) , (B7) 

Poc 2 



where the fractions / specify the energy partition fractions of Section 12721 Specifically, fa = Q/U m , /b = (M x + M z )/U m , f vx = 
K x /U m , f vz = K z /U m , and f xz = (fvxfvz) 1 ^ ■ We assume that protons dominate other ion species in the viscosity and therm al 
conductivity terms, and that ele ctrons dominate the electrical resistivity terms (see also iTulfl 984t IWhangll 1 9971; ICamposll 1999b . 

We use the Braginskii (1965) expressions for the transport c oefficients in the fully collisional limit. These coefficients depend 
on the proton and electron Coulomb collision timescales (e.g., Spitze3 [T962h . 

3 [^(WV 2 (B8) 



4 V 7T e 4 n p lnA 

3 /^7 M) 3/2 _„ 
Te = 4V^^TnA (B9) 

where we approximate the Coulomb logarithm in the coronal regime of temperatures and densities by 

In A = 23.2 + -In | — ^— )--lnf " e A . (BIO) 
2 V 10 6 K y 2 V10 6 cm-V V ' 

We also specify the magnitudes of the proton and electron gyrofrequencies, 

O p = f^,a = ^ (Bll) 
rtipC m e c 

and the dimensionless products 

x p = TpVLp , x e = T e Q, e . (B12) 

Thus, the proton viscosity coefficients are given by 

r)a = Q.96n p k B TpTp , rji = 2x p t] 2 (B13) 



, 1.24+2.23 , 

m=n D k Si T v T D —. . (B14) 

12 pappy ^ + 4.03*2 + 2.33 J 



The electrical conductivities are given by 
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FIG. 18. — Lineal' collisional damping rates forMHD waves. Color denotes the physical dissipation process: viscosity (black), electrical resistivity (red), and 
thermal conductivity (blue). Line style denotes the wave mode: Alfven (solid curves), fast mode (dashed curves), and slow mode (dotted curves). All quantities 
are plotted as base-10 logarithms of the rates, in units of s , for the standard model with parameters listed in Table[T] 



The thermal conductivities are given by 



«n = 



6.416^+1.837 
"4+14.794 + 3.7703 

3.906 n p k B T p T p 



(B16) 



(B17) 



n p k B T p T p 



2^ + 2.645 
4 + 2.704+0.677 



(B18) 



Finally, w e need to take acc ount of the transition from collisional to collisionless wave damping. In the low-density limit of 
the classical Braginskii (1965) expressions, some of the transport coefficients (e.g., 770, rj\, kii) become infini tely large as the 
mean time between collisions beco mes infinite. This "molasses limit" has been recognized to be unphvs ical ("see | Williamsll 1995); 
Cranmer & van Ballegooiien 2005). Thus, we derived a simplified version of the general expressions of Chang & Callen (1992) 
and lJi et all (12009b to describe what happens when collisions become infrequent. We computed the 7 dissipation rates as above, 
but then we multiplied them all by the following dimensionless factor C, where 



C- 



7"esc + 7~n 



and the macroscopic expansion timescale for waves is estimated as 



(u + V A )\dp /dr 



(B19) 



(B20) 



For strong collisions, C« 1 and the Bragi nskill (1 19651) expressions are valid. For weak collisions, C s» t bsc /t p -C 1 and the 
damping rates are quenched. 

Figure[l8]illustrates how the components of the wave damping rates vary with radial distance in the model of the fast solar wind 
described in Section [2] For the fast and slow mode waves, the viscous and conductive damping terms are of roughly comparable 
strength, but for the fast mode the conductive damping wins out at large distances (f3 3> 1). The viscous term is most important 
for the Alfven mode, but its magnitude remains small in comparison to the dominant terms for fast and slow mode damping. At 
the heights displayed here, Ohmic dissipation never appears t o be important in compari son to the other terms. This situation is 
reversed, however, lower down in the chromosphere (see, e.g..lKhodache nko et alJ l2004). The curves in Figure[l8]are shown for 
the general case of the transition to a collisionless plasma ( i.e., with all rates multiplied by C). For z < 0.1 Rq in the low corona, 
C~l and the general rates are identical to the unmodified BraginskiJ (1 19651) rates. For heights greater than z ~ 1^©, however, 
the the rates multiplied by C become about two orders of magnitude smaller than the unmodified classical rates. 
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C. ANALYTIC SOLUTIONS TO ADVECTION-DIFFUSION PROBLEMS IN LIMITED PARAMETER REGIMES 

C.l. Alfven Waves: Cascade and Source Terms 

Equation ( T38T > is a reduced one-dimensional v ersion of the full advection-diffusion equation for Alfvenic fluctuations. The 
Appendix 

ofE 

ranmer & van B allegooiien (2003) presented one method of solving this equation in the low-wavenumber, strong 
turbulence (xo "C 1) limit. Here we derive a more general case for arbitrary xo- Ignoring both wave damping and mode coupling, 
and assuming a steady state (i.e., db\/dt = 0), Equation (1381 can be simplified to 

7T = §A (C1) 

ox 

where here we define x = lnk± and we write the cascade rate as 

If - *i\ 
e=— Uj&a-OL— ± (C2) 



r A V dx 
a±_k l l s d , 2 



... yy±kl) (C3) 

where s = fj,±/a±. Note that the ratio s was called f3/j by ICranmer & van B allegooiienl (120031) and lLandi & Cranmerl (f2009). 
The second form for e given in Equation (IC3b helps to show that the power-law spectrum for b± in the inertial range (i.e., where 
§\ = and e is constant) should be independent of the value of s. In the limiting cases of strong (xo < 1) and weak (xo 3> 1) 
turbulence, we use Equation d36*l l to find that b± is proportional to k^ 3 and k^ 2 , respectively. 

In regions of wavenumber space where the source term is nonzero, e is not constant and the simple inertial-range scalings do 
not apply. If we assume that most of the fluctuation power is injected near i«io = ln&oj_, then it makes sense to use a compact 
Gaussian shape for the source term, 



Sa(x) = J° exp 



x—xo 

&0 



(C4) 



where the dimensionless width of the Gaussian is spe cified by op = 1 in our models. Th e constant £o is varied arbitrarily to 
produce the desired total fluctuation energy density U\ . ICranmer & van Ballegooijenl (120031) showed how the above form for the 
source function integrates to a cascade rate 

' x— xq s 



six) = — 



1+erf 



Finally, we define an auxiliary parameter q = b 2 ± k~^ and integrate Equation ( !C3t to obtain 

3/2 _ 3 f°° dn(l + xo)e(K) 



(C5) 



q 2a±J k± 0V2 K 2 + (3. s /2) ■ < C6) 

In practice, we integrate this equation numerically and use the definition of q to obtain b±(k±). In the energy containing range 
(x <C Xq), we see that s — > 0, and thus q is constant and b\ oc k s ± . The low-^j^ region of wavenumber space stands in contrast to 
the inertial range because here the shape of the fluctuation spectrum does depend on the value of s. 
In the MHD strong-turbulence inertial range (where e « so, <f> ~ 1, and xo 1), we obtain the standard solution b± ocv± oc 

k^ 3 . However, when k^p p increases past unity into the KAW dispersive range, we see from Equation (l35l ) that for /3 3> 1 there 
exists a sizable "dispersion range" in which <f> oc k\ and we obtain 

b± oc £ J 2/ ' 3 , oc k + ^ 3 . (C7) 

Converting these into the more commonly used one-dimensional spectra (see AppendixlAl. the MHD inertial range has e A oc k^ 3 , 

—1 /3 —1/3 

and the KAW inertial range has e A oc k ± for the magnetic fluc t uations and e a oc k , for the (elect ron) velocity fluctuations . 
These scalings have been described by, e.g., Biskamp et al. ( 1996), Cranm er & van B allegooiien (2003), and lHowes et al.l d2008). 
Note, however, that strong damping also begins to occur in the KAW regime, so the above power laws may not be evident in the 
final modeled spectra. 

C.2. Alfven Waves: Cascade and Dissipation Terms 

We include the effects of high-fe^ dissipation in the perpendicular cascade by assuming the damping acts only at wavenumbers 
significantly above those where the source term is dominant. Thus, we assume that §a = and we continue to ignore mode 
coupling. In this limiting case, the time-steady advection-diffusion equation becomes 

k ± ^=-2%bl. (C8) 
ok± 

We note that lHowes et alj (120081) found it is important to include KAW damping when solving for the steady-state wave power 
(and thus the proton/electron energy partitioning) at large values of k±. 

In Appendix IC.ll we showed that the properti es of the inert i al ran ge spectra should be independent of the value of s. In order to 
produce a closed-form solution, here we follow lHowes et ai1d200 8) and assume that s — >• oo (i.e., the cascade proceeds purely by 
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wavenumber advection). Thus, by assuming that a± = 0, we can use Equation (IC2b to write the b\ term on the right-hand side 
of Equation dC8b as 

(C9) 



if we also assume xo "C 1 in the high-Ai^ limit. To retain the most generality in cases when s is not infinitely large, we can use 
Equation (l4"TT i to replace fi± by fix* = /ij_ + 2a^/3. This may not be completely accurate for the cases of weakest advection (i.e., 
s w 0), but it is an improvement on ignoring the influence of the a± diffu sion term altogether . 

To simplify the modified transport equation, we take a further cue from lHowes et al.l (2008) and use the critical balance condi- 
tion ui ss k±v± to rewrite the equation as 

*t*._2_fl) (cio) 

where the ratio (7/0;) a is the output of the Vlasov-Maxwell dispersion analysis discussed in Section [5] For the KAW domain, 
this ratio is largely independent of kn and thus it can be treated as a function of k± only. When we solve this equation numerically, 
we start the integration at a low enough value of k± that the damping is negligibly small (i.e., at which e = £0). Thus, we integrate 
upwards in k± with 



e = £0 exp 





f dk± 











(Cll) 



Finally, the damped solution for e(k±) is used in Equation ( IC6b to obtain the damped power spectrum. 



C.3. Alfven Waves: Coupled Parallel and Perpendicular Transport 

The previous sections described the cascade as a function of k± and ignored the behavior of the full power spectrum 
Epjik\\,k±). To first order, the strong predicted anisotropy of MHD turbulence justifies this approach, but we are also con- 
cerned wit h the possible leakage of power to high values of k\\ and thus to high frequencies. Here we mainly follow the 
analysis of ICra nmer & van Balleg ooiienl (120031) . but we also include the possible effects of weak turbulence when \o 3> 1- 
Goldreich & Sridhaii (ri995l) wrote the full power spectrum as a separable function of two variables: kj_ and %, with 

V A b±ik±) 

'A 

and x = k\\V\/(k±b±). This definition allows the dimensionless function g(x) to be normalized to unity, 



E A (k h k ± )= , 3 - Six) (C12) 

K 1 



d X g(x) = 1 . (CO) 

but in practice we usually calculate the normalization for g(x) from the condition that the total power over all wavenumber space 
integrates properly to Ua- 

It has been known for some time that the dominant contribution to the integral in Equation ( |C13t should come from the region 

where \x\ < 1. For values of ku at which \x\ < 1, the iGoldreich & Sridharl (1 19951) solution for b± oc k^_ gives a dominant 

perpendicular dependence for the intertial range of E\ oc ^ i 10//3 . We also expect g(x) to grow neglig ibly small for |y | 3> 1, and 
thus for these large-fc|| regions of wavenumber space there should be very little Alfvenic wave power. ICho et al.l (|2002) found that 
numerical simulations of anisotropic MHD turbulence were consistent with g{x) being fit reason ably well with either a simple 
expon ential function (g ~ e~ x ) or a Castaing function (a convolution of multiple exponentials). ICranmer & van B allegooiien 
(2003) derived an analytic solution to a simplified version of Equation OTb . with 

W . «» f,^ (C1 4) 



3T(n- 0.5)^7? V 9a 
and 

3s 

n = l + — . (C15) 

-1/3 

These expressions are appropriate for the MHD inertial range where vj_ oc k ± , but we use them for the ent ire range of mod- 
eled wavenumbers. Th e above form for g(x) resembles a generalized Lorentzian, or kappa distribution (e.g., Vasvliunasl l 19681 : 
iPierrard & LazarlfeOlOl) that is Gaussian for small arguments and evolves to a power-law tail for large arguments. We can 
simplify the argument of the power-law term by using the values of the cascade parameters discussed in Section 13.11 with 
a±/ai\ w 18.6/(35+2). 

An example choice for the dimensionless advection-diffusion ratio (s = fi± /a± = 2) gives rise to a power-law exponent n = 5/2, 
and the large-£|| behavior of the Alfvenic power spectrum is E A o c fcjj" 5 '. Smaller values of s give shallower power-law slopes. In 
fact. ICranmer & van Ballegooiien (2003} and lLandi & Cranmerl (120091) found that if s could be maintained at small values of 
order 0.1-0.3, there would be sufficient high-fc[| power to heat protons and heavy ions in the corona via ion cyclotron resonance. 
In the opposite limit of pure advection (i.e., i^ooor a± -+ 0, with an and fi± remaining finite) Equation ( IC14l i becomes a 
Gaussian, 

gtxJocexpf-^x 2 ). (C16) 
V 3aii / 
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IChandrarJ d2008bl) also obtained a similar Gaussian solution for the parallel spectrum under the assumption of pure advection. 
Using the values of the cascade parameters discussed in Section lXTl we can set [i± w 1 .95 in the limit of s — > oo. Thus, the ratio 
[ix/atu ~ 6.2 and we can write g ~ e~ 2x in the pure-advection limit. 

In this paper we modify the analysis described above in one additional way. Instead of using the usual critical balance parameter 
X as the argument of g(x)> we instead use 

Xeff = JL- (C17) 



1 + Xo 

where xo is defined in Equation yT}. In the strong turbulence regime (xo •€ 1) this modification makes no difference. In the 
weak turbulence regime (xo 3> 1) this has the effect of extending the "filled" region of the spectrum (i.e., g(Xeff) ~ 1) up through 
all wavenumbers with kn < & || ■ 

Finally, we must adjust the highest frequency part of the spectrum to account self-consistently for the effects of ion cyclotron 
damping. Because the cyclotron resonance at high k\\ has a rapid onset (see Figure ITTT a)), we need only model its effects over a 
limited range of wavenumber space. We truncate the calculation of the spectrum at a maximum parallel wavenumber 

*»~ Va - 0J2 (CIS) 



n p 

at which Ita/O^I ~ L Abov e this wave number, we found that slowly-varying solutions to the Vlasov-Maxwell dispersion 
relation cease to exist (see also Stix 1992). Between 0.1£|| max and fc|| max , we include the time-steady effect of resonant damping 
by assuming that the local Alfvenic wave power is produced solely by the nonlinear coupling with the fast mode. If only the 
coupling and damping are present, the time-steady transport equation simplifies to 

OE a Ep — E a 

-5T ~ — " 2 7a£a = , (C19) 

at T AF 

which can be solved analytically for However, we note that we already have a time-steady solution for £a in the presence of 
nonlinear coupling, but it does not take into account the ion cyclotron damping. Equation ( 154-b gives that solution, which we now 
call Eqa- Thus, we insert it in place of Ep in Equation (IC19b above, since that is the solution toward which the coupling will drive 
the system in the absence of damping. We then use the analytic solution 

r ' - ti)x (C20) 



1+27ATAF 

to account for ion cyclotron dissipation at high k\\ . This solution gives rise to a significant reduction in the power spectrum when 
the damping rate 7a exceeds the rate at which power is supplied from the nonlinear coupling. 

C.4. Fast-Mode Waves: Cascade and Source Terms 

This section is conceptually similar to Appendix lC.il in that we ignore both damping and mode coupling, assume time-steady 
conditions, and model the spectral transport of fast-mode waves as a balance between diffusive cascade and the outer-scale source 
term. In that limiting case, Equation d32l becomes 

M = k2sF (C21) 



where the cascade rate is defined here as 

4-napk s sin0 dEp 
£ = — E *l*- (C22) 

When Sp = 0, the cascade rate is constant along radial rays of constant 9, and thus Ep oc kT 1 ! 2 . We follow Chandran (20051) and 
others by assuming a Gaussian shape for the fast-mode source term, but we also add a corresponding sin 9 angle dependence, 
with 

k 



Sp(k) = So exp 



koF 



sin9 . (C23) 



We eventually set 5o at the level required to maintain the spectrum at the known total energy density Up. Our choice to constrain 
kop to be equal to k^± = 1/ Aj_ is discussed in Section [3Tl 

Because both the left and right sides of Equation (lC211 i depend identically on sin 6, the resulting time-steady solution for 
Ep(k) becomes independent of 9. This outcom e was motivated by simulation results that sh ow that the fast-mode power spec- 
trum is largely isotropic in wavenumber space (Cho & Lazarian 2003^ ISvidzinski et al.ll2009l). It is also pos sible that additional 
isotropization of the fast-mode spectrum can come from couplings with slow-mode waves (Chandran 2008a J) or from multi -scale 
"wandering" of the magnetic field that gives rise to a continuously varying direction for 9 = (IShalchi & Kourakiil2007l) . In 
future work, our method of artificially forcing isotropy via the source term should be replaced with a more realistic description. 

In any case, we cancel out both instances of sin 9 and integrate Equation ( IC2U to obtain 

e(k) = k 3 pSo | ^Icrf.v-^- ) (C24) 




28 



CRANMER AND VAN BALLEGOOIJEN 



where here x = k/k^p. In the limit x ^> 1, the term in parentheses above approaches a constant value of \/7r/4. In the limit x <C 1, 
the term in parentheses is approximately equal to x 3 /3. Finally, we integrate the definition of e to get the time-steady spectrum, 

E l V A f°° e(K) 



2 47rap Jk K 

which we evaluate numerically using Equation dC24b for the cascade rate in the integrand. In the energy containing range 
(k <C &of), this yields Ep oc k~ 2 , and thus v# oc A;" 1 " 1 / 2 . 

C.5. Fast-Mode Waves: Cascade and Dissipation Terms 

If we consider high values of k above those affected by the outer-scale source term, we can solve for the transition from the 
inertial range to the dissipation range in the fast-mode power spectrum. An analytic solution becomes possible if we rewrite the 
fast-mode spectral transport time rp as a function of k and 9 only. This can be done by using the time-steady inertial range scaling 
for v*, with 

v* = 4 (~) , r F = 2 . ^ A • (C26) 
\ko J VySmoy/kkQ 

Note that Equation (8) of ISuzuki et al.l (120071) gave the same result for the fast-mode cascade timescale (but without the sin# 
term). The normalization wavenumber ko is defined arbitrarily here; it needs to be set well below the regime of strong damping, 
but well above the outer-scale wavenumber A:of so that we can justify ignoring the source term. 

The above approximation gives Dp oc k 5 ' 2 sin#. It is also straightforward to model the fast-mode damping rate as being 
proportional to a constant power of the wavenumber, and thus we assume 7f = ^{k/kof. Note that the normalizing constant 70 
may depend on the angle 9 as well. The time-steady version of Equation d32| i becomes 

avvlklJ 1 sin 9 d f Qll dEp \ ( k\~ 

- L i^^rv)= 2 ^u) £ - (c2?) 

Defining the auxiliary variable 



7/2 

v = = m (C28) 



2 fk 



helps to greatly simplify the differential equation. Thus, 

d 2 E F c 7 £f 



dx 2 ^(2z+13)/7 



(C29) 



where 



2 7 oV A /2\ (2z+13)/7 



afV^kosm9 



is a constant that is essentially the ratio of the damping rate to the cascade rate at the normalization wavenumber k(,. 

For z > 1, Equation ( |C29t is solved analytically with two linearly independent terms proportional to the two types of mod- 
ified Bessel function (/„ and K n ). Knowing that the only physically realistic solution is one that decreases monotonically with 
increasing k (or with decreasing x), we then use only one of those terms, which is given by 



E F (k) oc k~ V4 K ( 



2C 



2 \ 1/(20 fk\ 7/(40 



1 wl \koJ 



(C31) 



where £ = 7/ (2z- 1). Note that the transit-time damping rate of Equation (l58l gives z = 1 and C = 7. In the limiting case that the 
modified Bessel function of the second kind has a small argument, we have K^(x) ~ x~^ and thus E F oc k~ 7 / 2 , independent of the 
value of £. This is the proper inertial-range solution in the case of either low wavenumber (k -C ko) or weak damping (c~ <C 1). 
The opposite case of a large argument gives expone ntially steep dissipation in the limit of large k and/or large 70. This kind of 
solution was also derived by Hunana & Zankj d2010l) . 

Another useful special case for the damping exponent is z = 1/2. For this value of the exponent, Equation (IC29b is solved with 
two linearly independent power-law terms. As above, we keep only the solution that does not diverge as x —> (i.e., as k —> 00), 
and the time-steady spectrum is given by 

Erik) oc r 7(I+ V^;)/ 4 . (C32) 

The weak-damping limit of c 7 <C 1 gives the proper inertial-range solution E? oc fe" 7 ' 2 , but the presence of a nonzero value of c 7 
makes the spectrum steeper. This is one (possibly rare) case in which a physically motivated source of damping gives rise to a 
power-law "dissipation range." 
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